System and method for image segmentation in generating computer models of a joint to undergo arthroplasty

ABSTRACT

A custom arthroplasty guide and a method of manufacturing such a guide are disclosed herein. The method of manufacturing the custom arthroplasty guide includes: a) generating medical imaging slices of the portion of the patient bone; b) identifying landmarks on bone boundaries in the medical imaging slices; c) providing model data including image data associated with a bone other than the patient bone; d) adjusting the model data to match the landmarks; e) using the adjusted model data to generate a three dimensional computer model of the portion of the patient bone; f) using the three dimensional computer model to generate design data associated with the custom arthroplasty guide; and g) using the design data in manufacturing the custom arthroplasty guide.

CROSS REFERENCE TO RELATED APPLICATION

The present application is a continuation application of U.S.application Ser. No. 15/581,974 filed Apr. 28, 2017, which applicationis a continuation of U.S. application Ser. No. 14/946,106 filed Nov. 19,2015, now U.S. Pat. No. 9,687,259, which application is a continuationof U.S. application Ser. No. 13/731,697 filed Dec. 31, 2012, now U.S.Pat. No. 9,208,263, which application is a continuation of U.S.application Ser. No. 13/374,960 filed Jan. 25, 2012, now U.S. Pat. No.8,532,361, which application is a continuation of U.S. patentapplication Ser. No. 13/066,568, filed Apr. 18, 2011, now U.S. Pat. No.8,160,345, which application is a continuation-in-part application ofU.S. patent application Ser. No. 12/386,105 filed Apr. 14, 2009, nowU.S. Pat. No. 8,311,306. U.S. application Ser. No. 12/386,105 claims thebenefit under 35 U.S.C. § 119(e) of U.S. Provisional Patent ApplicationNo. 61/126,102, entitled “System and Method For Image Segmentation inGenerating Computer Models of a Joint to Undergo Arthroplasty” filed onApr. 30, 2008. Each of these applications is hereby incorporated byreference herein in their entireties.

FIELD OF THE INVENTION

The present invention relates to image segmentation. More specifically,the present invention relates to image segmentation in generatingcomputer models of a joint to undergo arthroplasty, wherein the computermodels may be used in the design and manufacture of arthroplasty jigs.

BACKGROUND OF THE INVENTION

Over time and through repeated use, bones and joints can become damagedor worn. For example, repetitive strain on bones and joints (e.g.,through athletic activity), traumatic events, and certain diseases(e.g., arthritis) can cause cartilage in joint areas, which normallyprovides a cushioning effect, to wear down. Cartilage wearing down canresult in fluid accumulating in the joint areas, pain, stiffness, anddecreased mobility.

Arthroplasty procedures can be used to repair damaged joints. During atypical arthroplasty procedure, an arthritic or otherwise dysfunctionaljoint can be remodeled or realigned, or an implant can be implanted intothe damaged region. Arthroplasty procedures may take place in any of anumber of different regions of the body, such as a knee, a hip, ashoulder, or an elbow.

One type of arthroplasty procedure is a total knee arthroplasty (“TKA”),in which a damaged knee joint is replaced with prosthetic implants. Theknee joint may have been damaged by, for example, arthritis (e.g.,severe osteoarthritis or degenerative arthritis), trauma, or a raredestructive joint disease. During a TKA procedure, a damaged portion inthe distal region of the femur may be removed and replaced with a metalshell, and a damaged portion in the proximal region of the tibia may beremoved and replaced with a channeled piece of plastic having a metalstem. In some TKA procedures, a plastic button may also be added underthe surface of the patella, depending on the condition of the patella.

Implants that are implanted into a damaged region may provide supportand structure to the damaged region, and may help to restore the damagedregion, thereby enhancing its functionality. Prior to implantation of animplant in a damaged region, the damaged region may be prepared toreceive the implant. For example, in a knee arthroplasty procedure, oneor more of the bones in the knee area, such as the femur and/or thetibia, may be treated (e.g., cut, drilled, reamed, and/or resurfaced) toprovide one or more surfaces that can align with the implant and therebyaccommodate the implant.

Accuracy in implant alignment is an important factor to the success of aTKA procedure. A one- to two-millimeter translational misalignment, or aone- to two-degree rotational misalignment, may result in imbalancedligaments, and may thereby significantly affect the outcome of the TKAprocedure. For example, implant misalignment may result in intolerablepost-surgery pain, and also may prevent the patient from having full legextension and stable leg flexion.

To achieve accurate implant alignment, prior to treating (e.g., cutting,drilling, reaming, and/or resurfacing) any regions of a bone, it isimportant to correctly determine the location at which the treatmentwill take place and how the treatment will be oriented. In some methods,an arthroplasty jig may be used to accurately position and orient afinishing instrument, such as a cutting, drilling, reaming, orresurfacing instrument on the regions of the bone. The arthroplasty jigmay, for example, include one or more apertures and/or slots that areconfigured to accept such an instrument.

A system and method has been developed for producing customizedarthroplasty jigs configured to allow a surgeon to accurately andquickly perform an arthroplasty procedure that restores thepre-deterioration alignment of the joint, thereby improving the successrate of such procedures. Specifically, the customized arthroplasty jigsare indexed such that they matingly receive the regions of the bone tobe subjected to a treatment (e.g., cutting, drilling, reaming, and/orresurfacing). The customized arthroplasty jigs are also indexed toprovide the proper location and orientation of the treatment relative tothe regions of the bone. The indexing aspect of the customizedarthroplasty jigs allows the treatment of the bone regions to be donequickly and with a high degree of accuracy that will allow the implantsto restore the patient's joint to a generally pre-deteriorated state.However, the system and method for generating the customized jigs oftenrelies on a human to “eyeball” bone models on a computer screen todetermine configurations needed for the generation of the customizedjigs. This “eyeballing” or manual manipulation of the bone modes on thecomputer screen is inefficient and unnecessarily raises the time,manpower and costs associated with producing the customized arthroplastyjigs. Furthermore, a less manual approach may improve the accuracy ofthe resulting jigs.

There is a need in the art for a system and method for reducing thelabor associated with generating customized arthroplasty jigs. There isalso a need in the art for a system and method for increasing theaccuracy of customized arthroplasty jigs.

SUMMARY

Systems and methods for image segmentation in generating computer modelsof a joint to undergo arthroplasty are disclosed. Some embodiments mayinclude a method of partitioning an image of a bone into a plurality ofregions, where the method of partitioning may include obtaining aplurality of volumetric image slices of the bone, generating a pluralityof spline curves associated with the bone, verifying that at least oneof the plurality of spline curves follow a surface of the bone, andcreating a three dimensional (3D) mesh representation based upon the atleast one of the plurality of spline curves.

Other embodiments may include a method of generating a representation ofa model bone, where the method of generating the representation mayinclude obtaining an image scan of the representation as a plurality ofslices, segmenting each slice in the plurality into one or moresegmentation curves, generating a mesh of the representation, adjustingeach slice in the plurality to include areas where the contact area ofthe bone is stable between successive image scans, and generating anchorsegmentation such that the anchor segmentation follows a boundary of therepresentation of the model bone.

Other embodiments may include a method of segmenting a target bone usinga representation of a model bone, where the method of segmenting thetarget bone may include registering a segmented form of therepresentation to an image scan of the target bone, refining theregistration of the segmented form of the representation near a boundaryof the target bone, generating a mesh from the segmented form of therepresentation, and generating a plurality of spline curves thatapproximate the intersection of the mesh and one or more slices from theimage scan of the target bone.

Other embodiments may include a method of mapping a representation of amodel bone into an image scan of a target bone, where the method ofmapping may include registering a generated portion of therepresentation into the image scan of the target bone using atranslational transformation, registering the generated portion of therepresentation into the image scan of the target bone using a similaritytransformation, registering a boundary portion of the representationinto the image scan of the target bone using an affine transformation,and registering the boundary portion of the representation into theimage scan of the target bone using a spline transformation.

Other embodiments may include a method for determining a degree ofcorrespondence between an image of a target bone and a representation ofa model bone, where the method of determining correspondence may includeselecting a plurality of sample points in the representation of themodel bone to be registered, partitioning the plurality of sample pointsinto a plurality of groups, sampling the image of the target bone,determining a correlation of voxel intensities between the image of thetarget bone and the representation of the model bone for each group inthe plurality, and averaging the correlation determined for each groupin the plurality.

Other embodiments may include a method for refining registration of arepresentation of a model bone to a target bone, where the method ofrefining may include transforming an anchor segmentation mesh,generating a plurality of random points around the transformed anchorsegmentation mesh, determining if each point in the plurality liesinside one or more of the following meshes: InDark-OutLight,InLight-OutDark, or Dark-In-Light, determining whether one or more ofthe plurality of points lie within a threshold distance of the surfaceof the transformed anchor segmentation mesh, and adding each point inthe plurality as a dark point or light point depending upon whether thepoint lies within the InDark-OutLight, InLight-OutDark, or Dark-In-Lightmeshes.

Still other embodiments may include a method for generating splinecurves outlining the surface of a feature of interest of a target bone,where the method of generating spline curves may include intersecting a3D mesh model of the feature surface with one or more slices of targetdata (the intersection defining a polyline curve), parameterizing thepolyline curve as a function of length and tangent variation,calculating a weighted sum of the length and tangent parameterizations,and sampling the polyline using the results of the act of calculating.

A custom arthroplasty guide and a method of manufacturing such a guideare disclosed herein. The guide manufactured includes a mating regionconfigured to matingly receive a portion of a patient bone associatedwith an arthroplasty procedure for which the custom arthroplasty guideis to be employed. The mating region includes a surface contour that isgenerally a negative of a surface contour of the portion of the patientbone. The surface contour of the mating region is configured to matewith the surface contour of the portion of the patient bone in agenerally matching or interdigitating manner when the portion of thepatient bone is matingly received by the mating region. The method ofmanufacturing the custom arthroplasty guide includes: a) generatingmedical imaging slices of the portion of the patient bone; b)identifying landmarks on bone boundaries in the medical imaging slices;c) providing model data including image data associated with a boneother than the patient bone; d) adjusting the model data to match thelandmarks; e) using the adjusted model data to generate a threedimensional computer model of the portion of the patient bone; f) usingthe three dimensional computer model to generate design data associatedwith the custom arthroplasty guide; and g) using the design data inmanufacturing the custom arthroplasty guide.

In one embodiment of the method of manufacturing the custom arthroplastyguide, operation a) includes generating at least one of MRI slices or CTslices of the portion of the patient bone.

In one embodiment of the method of manufacturing the custom arthroplastyguide, operation b) includes placing landmark points on the boneboundaries in the medical imaging slices. For example, placing thelandmark points may include a user at a user interface employing atleast one of a mouse, keyboard, pen-and-tablet system, touch screensystem, or spatial input device to place landmark points. The boneboundaries may include lines representative in the medical imagingslices of cortical bone boundaries.

In one embodiment of the method of manufacturing the custom arthroplastyguide, the image data associated with a bone other than the patient bonemay include a golden bone model. In one embodiment of the method ofmanufacturing the custom arthroplasty guide, the image data of operationc) may include a model of a bone that is the same bone type as thepatient bone, but representative of what is considered to be generallynormal with respect to the bone type for at least one of size,condition, shape, weight, age, height, race, gender or diagnosed diseasecondition.

In one embodiment of the method of manufacturing the custom arthroplastyguide, operation d) includes deforming a mesh to match the landmarks,wherein the mesh is associated with the model data. For example,deforming the mesh to match the landmarks may include registering themesh to the landmarks using at least one of translation transforms,similarity transforms or affine transforms.

Deforming the mesh to match the landmarks may further include detectingimage edges near the mesh. For example, detecting image edges near themesh may include computing a signed distance image for the mesh.Detecting the image edges near the mesh may further include computing agradient of the signed distance image. Detecting the image edges nearthe mesh may yet further include computing a gradient of at least one ofthe medical imaging slices. Detecting the image edges near the mesh maystill further include computing an edges image by correcting thegradient of at least one of the medical imaging slices with the gradientof the signed distance image.

In one embodiment of the method of manufacturing the custom arthroplastyguide, deforming the mesh to match the landmarks may yet further includeregistering simultaneously the mesh to the landmarks and to the imageedges using B-spline deformable transforms.

In one embodiment of the method of manufacturing the custom arthroplastyguide, operation d) may further include generating segmentation curvesfrom the deformed mesh. Operation d) may yet further includeapproximating the segmentation curves with spline curves. Additionally,operation d) may still further include modifying the spline curves tomatch the landmarks.

In one embodiment of the method of manufacturing the custom arthroplastyguide, modifying the spline curves to match the landmarks may includeidentifying a landmark from the landmarks and computing a distance ofthe spline curve to the indentified landmark. Modifying the splinecurves to match the landmarks may further include identifying a closestarc of the spline curve to the identified landmark and modifying theclosest arc to include the identified landmark, resulting in a modifiedspline curve. Modifying the spline curves to match the landmarks may yetfurther include computing distances from other landmarks to the modifiedspline curve and, if the distances are worse, disregarding themodification to the closest arc that resulted in the modified splinecurve.

If the modification is disregarded, then modifying the spline curves tomatch the landmarks may still further include identifying a point in theclosest arc that is closest to the identified landmark and inserting aspline vertex at the identified point. Modifying the spline curves tomatch the landmarks may still further include moving the inserted vertexto the identified landmark.

In one embodiment of the method of manufacturing the custom arthroplastyguide, in operation e) the three dimensional computer model of theportion of the patient bone may include: a bone only model thatrepresents bone surfaces but no cartilage surfaces; or a restored bonemodel that represents bone surfaces but no cartilage surfaces.

In one embodiment of the method of manufacturing the custom arthroplastyguide, in operation e) the three dimensional computer model of theportion of the patient bone includes an arthritic model that representsboth bone and cartilage surfaces.

In one embodiment of the method of manufacturing the custom arthroplastyguide, in operation f) the data associated with the custom arthroplastyguide may include at least one of data associated with the design of thesurface contour of the mating region or a resection plane associatedwith a resection guide surface of the custom arthroplasty guide.

In one embodiment of the method of manufacturing the custom arthroplastyguide, in operation f) the data associated with the custom arthroplastyguide includes manufacturing instructions for manufacturing the customarthroplasty guide.

In one embodiment of the method of manufacturing the custom arthroplastyguide, operation g) may include using the design data to guide a CNCmachine or SLA machine (or other analogous additive or subtractivemanufacturing technologies) in manufacturing the custom arthroplastyguide.

Yet another embodiment of manufacturing the above-described customarthroplasty guide is disclosed herein. Specifically, the method ofmanufacturing the custom arthroplasty guide includes: a) generatingmedical imaging data of the portion of the patient bone; b) identifyinglandmarks associated with bone boundaries in the medical imaging data;c) providing a geometrical representation of an exemplary bone that isnot the patient bone but is the same bone type as the patient bone andis representative of what is considered to be generally normal withrespect to the bone type for at least one of size, condition, shape,weight, age, height, race, gender or diagnosed disease condition; d)adjusting the geometrical representation of the exemplary bone to matchthe landmarks; e) using the adjusted geometrical representation of theexemplary bone to generate a three dimensional computer model of theportion of the patient bone; f) using the three dimensional computermodel to generate design data associated with the custom arthroplastyguide; and g) using the design data in manufacturing the customarthroplasty guide. In one version of the embodiment, the geometricalrepresentation of an exemplary bone includes a golden femur mesh or agolden tibia mesh.

Yet still another embodiment of manufacturing the above-described customarthroplasty guide is disclosed herein. Specifically, the method ofmanufacturing the custom arthroplasty guide includes: a) generatingmedical imaging slices of the portion of the patient bone; b)identifying landmarks on cortical bone boundaries in the medical imagingdata; c) providing golden bone image data that is not of the patientbone but is of the same bone type as the patient bone and isrepresentative of what is considered to be generally typical withrespect to a characteristic of bones such as the patient bone; d)segmenting the golden bone image data and generating a golden bone meshfrom the segmented golden bone image data; e) deforming the golden bonemesh to match the landmarks; f) generating at least one of segmentationcurves or spline curves from the deformed golden bone mesh; g)(optional) modify the at least one of segmentation curves or splinecurves to match the landmarks; h) using the modified at least one ofsegmentation curves or spline curves to generate a three dimensionalcomputer model generally representative of the portion of the patientbone; i) using the three dimensional computer model to generate designdata associated with the custom arthroplasty guide; and j) using thedesign data in manufacturing the custom arthroplasty guide. In someversions of the embodiment with respect to operation c), thecharacteristic includes at least one of size, condition, shape, weight,age, height, race, gender or diagnosed disease condition.

While multiple embodiments are disclosed, still other embodiments of thepresent invention will become apparent to those skilled in the art fromthe following detailed description, which shows and describesillustrative embodiments of the invention. As will be realized, theinvention is capable of modifications in various aspects, all withoutdeparting from the spirit and scope of the present invention.Accordingly, the drawings and detailed description are to be regarded asillustrative in nature and not restrictive.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a schematic diagram of a system for employing the automatedjig production method disclosed herein.

FIGS. 1B-1E are flow chart diagrams outlining the jig production methoddisclosed herein.

FIGS. 1F and 1G are, respectively, bottom and top perspective views ofan example customized arthroplasty femur jig.

FIGS. 1H and 1I are, respectively, bottom and top perspective views ofan example customized arthroplasty tibia jig.

FIG. 2A is a sagittal plane image slice depicting a femur and tibia andneighboring tissue regions with similar image intensity.

FIG. 2B is a sagittal plane image slice depicting a region extendinginto the slice from an adjacent image slice.

FIG. 2C is a sagittal plane image slice depicting a region of a femurthat is approximately tangent to the image slice.

FIG. 3A is a sagittal plane image slice depicting an intensity gradientacross the slice.

FIG. 3B is a sagittal plane image slice depicting another intensitygradient across the slice.

FIG. 3C is a sagittal plane image slice depicting another intensitygradient across the slice.

FIG. 4A depicts a sagittal plane image slice with a high noise level.

FIG. 4B depicts a sagittal plane image slice with a low noise level.

FIG. 5 is a sagittal plane image slice of a femur and tibia depictingregions where good definition may be needed during automaticsegmentation of the femur and tibia.

FIG. 6 depicts a flowchart illustrating one method for automaticsegmentation of an image modality scan of a patient's knee joint.

FIG. 7A is a sagittal plane image slice of a segmented femur.

FIG. 7B is a sagittal plane image slice of a segmented femur and tibia.

FIG. 7C is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7D is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7E is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7F is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7G is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7H is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7I is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7J is another sagittal plane image slice of a segmented femur andtibia.

FIG. 7K is another sagittal plane image slice of a segmented femur andtibia.

FIG. 8 is a sagittal plane image slice depicting automatically generatedslice curves of a femur and a tibia.

FIG. 9 depicts a 3D mesh geometry of a femur.

FIG. 10 depicts a 3D mesh geometry of a tibia.

FIG. 11 depicts a flowchart illustrating one method for generating agolden template.

FIG. 12A is a sagittal plane image slice depicting a contour curveoutlining a golden tibia region, a contour curve outlining a grown tibiaregion and a contour curve outlining a boundary golden tibia region.

FIG. 12B is a sagittal plane image slice depicting a contour curveoutlining a golden femur region, a contour curve outlining a grown femurregion and a contour curve outlining a boundary golden femur region.

FIG. 13A depicts a golden tibia 3D mesh.

FIG. 13B depicts a golden femur 3D mesh.

FIG. 14A is a sagittal plane image slice depicting anchor segmentationregions of a tibia.

FIG. 14B is a sagittal plane image slice depicting anchor segmentationregions of a femur.

FIG. 15A is a 3D mesh geometry depicting the anchor segmentation mesh,the InDark-OutLight anchor mesh, the InLight-OutDark anchor mesh, andthe Dark-In-Light anchor mesh of a tibia.

FIG. 15B is a 3D mesh geometry depicting the anchor segmentation mesh,the InDark-OutLight anchor mesh and the InLight-OutDark anchor mesh of afemur.

FIG. 16 depicts a flowchart illustrating one method for performingautomatic segmentation of scan data using golden template registration.

FIG. 17 depicts a flowchart illustrating one method for mapping thesegmented golden femur template regions into the target scan data usingimage registration techniques.

FIG. 18 depicts a registration framework that may be employed by oneembodiment.

FIG. 19 depicts a flowchart illustrating one method for mapping thesegmented golden tibia template regions into the target scan data usingimage registration techniques.

FIG. 20 depicts a flowchart illustrating one method for computing ametric for the registration framework of FIG. 18.

FIG. 21 depicts a flowchart illustrating one method for refining theregistration results using anchor segmentation and anchor regions.

FIG. 22 depicts a set of randomly generated light sample points and darksample points of a tibia.

FIG. 23 depicts a flowchart illustrating one method for generatingspline curves to outline features of interest in each target MRI slice.

FIG. 24 depicts a polyline curve with n vertices.

FIG. 25 depicts a flowchart illustrating one method for adjustingsegments.

FIG. 26 is a sagittal plane image slice depicting a contour curve withcontrol points outlining a femur with superimposed contour curves of thefemur from adjacent image slices.

FIG. 27 depicts a 3D slice visualization of a femur showing the voxelsinside of the spline curves.

FIG. 28 is a diagram depicting types of data employed in the imagesegmentation algorithm that uses landmarks.

FIG. 29 is a flowchart illustrating the overall process for generating agolden femur model of FIG. 28.

FIG. 30 is an image slice of the representative femur to be used togenerate a golden femur mesh.

FIG. 31A is an isometric view of a closed golden femur mesh.

FIG. 31B is an isometric view of an open golden femur mesh created fromthe closed golden femur mesh of FIG. 31A.

FIG. 31C is the open femur mesh of FIG. 31B with regions of a differentprecision indicated.

FIGS. 32A-32B are isometric views of an open golden tibia mesh withregions of a different precision indicated.

FIG. 33 is a flow chart illustrating an alternative method of segmentingan image slice, the alternative method employing landmarks.

FIG. 34 is a flow chart illustrating the process involved in operation“position landmarks” of the flow chart of FIG. 33.

FIGS. 35A-35H are a series of sagittal image slices wherein landmarkshave been placed according the process of FIG. 34.

FIG. 36 is a flowchart illustrating the process of segmenting the targetimages that were provided with landmarks in operation “positionlandmarks” of the flow chart of FIG. 33.

FIG. 37 is a flowchart illustrating the process of operation “DeformGolden Femur Mesh” of FIG. 36, the process including mapping the goldenfemur mesh into the target scan using registration techniques.

FIG. 38A is a flowchart illustrating the process of operation “DetectAppropriate Image Edges” of FIG. 37.

FIG. 38B is an image slice with a contour line representing theapproximate segmentation mesh surface found in operation 770 c of FIG.37, the vectors showing the gradient of the signed distance for thecontour.

FIG. 38C is an enlarged view of the area in FIG. 38B enclosed by thesquare 1002, the vectors showing the computed gradient of the targetimage.

FIG. 38D is the same view as FIG. 38C, except the vectors of FIGS. 38Band 38C are superimposed.

FIG. 39 is a flowchart illustrating the process of operation “ModifySplines” of FIG. 36.

FIG. 40 is an image slice with a spline being modified according to theoperations of the flow chart of FIG. 39.

DETAILED DESCRIPTION

Disclosed herein are customized arthroplasty jigs 2 and systems 4 for,and methods of, producing such jigs 2. The jigs 2 are customized to fitspecific bone surfaces of specific patients. Depending on the embodimentand to a greater or lesser extent, the jigs 2 are automatically plannedand generated and may be similar to those disclosed in these three U.S.Patent Applications: U.S. patent application Ser. No. 11/656,323 to Parket al., titled “Arthroplasty Devices and Related Methods” and filed Jan.19, 2007; U.S. patent application Ser. No. 10/146,862 to Park et al.,titled “Improved Total Joint Arthroplasty System” and filed May 15,2002; and U.S. patent Ser. No. 11/642,385 to Park et al., titled“Arthroplasty Devices and Related Methods” and filed Dec. 19, 2006. Thedisclosures of these three U.S. Patent Applications are incorporated byreference in their entireties into this Detailed Description.

a. Overview of System and Method for Manufacturing CustomizedArthroplasty Cutting Jigs

For an overview discussion of the systems 4 for, and methods of,producing the customized arthroplasty jigs 2, reference is made to FIGS.1A-1E. FIG. 1A is a schematic diagram of a system 4 for employing theautomated jig production method disclosed herein. FIGS. 1B-1E are flowchart diagrams outlining the jig production method disclosed herein. Thefollowing overview discussion can be broken down into three sections.

The first section, which is discussed with respect to FIG. 1A and[blocks 100-125] of FIGS. 1B-1E, pertains to an example method ofdetermining, in a three-dimensional (“3D”) computer model environment,saw cut and drill hole locations 30, 32 relative to 3D computer modelsthat are termed restored bone models 28 (also referenced as “planningmodels” throughout this submission.) In some embodiments, the resulting“saw cut and drill hole data” 44 is referenced to the restored bonemodels 28 to provide saw cuts and drill holes that will allowarthroplasty implants to restore the patient's joint to itspre-degenerated state. In other words, in some embodiments, thepatient's joint may be restored to its natural alignment, whethervalgus, varus or neutral.

While many of the embodiments disclosed herein are discussed withrespect to allowing the arthroplasty implants to restore the patient'sjoint to its pre-degenerated or natural alignment state, many of theconcepts disclosed herein may be applied to embodiments wherein thearthroplasty implants restore the patient's joint to a zero mechanicalaxis alignment such that the patient's knee joint ends up being neutral,regardless of whether the patient's predegenerated condition was varus,valgus or neutral. For example, as disclosed in U.S. patent applicationSer. No. 12/760,388 to Park et al., titled “Preoperatively Planning AnArthroplasty Procedure And Generating A Corresponding Patient SpecificArthroplasty Resection Guide”, filed Apr. 14, 2010, and incorporated byreference into this Detailed Description in its entirety, the system 4for producing the customized arthroplasty jigs 2 may be such that thesystem initially generates the preoperative planning (“POP”) associatedwith the jig in the context of the POP resulting in the patient's kneebeing restored to its natural alignment. Such a natural alignment POP isprovided to the physician, and the physician determines if the POPshould result in (1) natural alignment, (2) mechanical alignment, or (3)something between (1) and (2). The POP is then adjusted according to thephysician's determination, the resulting jig 2 being configured suchthat the arthroplasty implants will restore the patient's joint to (1),(2) or (3), depending on whether the physician elected (1), (2) or (3),respectively. Accordingly, this disclosure should not be limited tomethods resulting in natural alignment only, but should, whereappropriate, be considered as applicable to methods resulting in naturalalignment, zero mechanical axis alignment or an alignment somewherebetween natural and zero mechanical axis alignment.

The second section, which is discussed with respect to FIG. 1A and[blocks 100-105 and 130-145] of FIGS. 1B-1E, pertains to an examplemethod of importing into 3D computer generated jig models 38 3D computergenerated surface models 40 of arthroplasty target areas 42 of 3Dcomputer generated arthritic models 36 of the patient's joint bones. Theresulting “jig data” 46 is used to produce a jig customized to matinglyreceive the arthroplasty target areas of the respective bones of thepatient's joint.

The third section, which is discussed with respect to FIG. 1A and[blocks 150-165] of FIG. 1E, pertains to a method of combining orintegrating the “saw cut and drill hole data” 44 with the “jig data” 46to result in “integrated jig data” 48. The “integrated jig data” 48 isprovided to the CNC machine 10 for the production of customizedarthroplasty jigs 2 from jig blanks 50 provided to the CNC machine 10.The resulting customized arthroplasty jigs 2 include saw cut slots anddrill holes positioned in the jigs 2 such that when the jigs 2 matinglyreceive the arthroplasty target areas of the patient's bones, the cutslots and drill holes facilitate preparing the arthroplasty target areasin a manner that allows the arthroplasty joint implants to generallyrestore the patient's joint line to its pre-degenerated state.

As shown in FIG. 1A, the system 4 includes one or more computers 6having a CPU 7, a monitor or screen 9 and an operator interface controls11. The computer 6 is linked to a medical imaging system 8, such as a CTor MRI machine 8, and a computer controlled machining system 10, such asa CNC milling machine 10.

In another embodiment, rather than using a single computer for the wholeprocess, multiple computers can perform separate steps of the overallprocess, with each respective step managed by a respective technicianskilled in that particular aspect of the overall process. The datagenerated in one process step on one computer can be then transferredfor the next process step to another computer, for instance via anetwork connection.

As indicated in FIG. 1A, a patient 12 has a joint 14 (e.g., a knee,elbow, ankle, wrist, hip, shoulder, skull/vertebrae orvertebrae/vertebrae interface, etc.) to be replaced. The patient 12 hasthe joint 14 scanned in the imaging machine 8. The imaging machine 8makes a plurality of scans of the joint 14, wherein each scan pertainsto a thin slice of the joint 14.

As can be understood from FIG. 1B, the plurality of scans is used togenerate a plurality of two-dimensional (“2D”) images 16 of the joint 14[block 100]. Where, for example, the joint 14 is a knee 14, the 2Dimages will be of the femur 18 and tibia 20. The imaging may beperformed via CT or MRI. In one embodiment employing MRI, the imagingprocess may be as disclosed in U.S. patent application Ser. No.11/946,002 to Park, which is entitled “Generating MRI Images Usable ForThe Creation Of 3D Bone Models Employed To Make Customized ArthroplastyJigs,” was filed Nov. 27, 2007 and is incorporated by reference in itsentirety into this Detailed Description.

As can be understood from FIG. 1A, the 2D images are sent to thecomputer 6 for creating computer generated 3D models. As indicated inFIG. 1B, in one embodiment, point P is identified in the 2D images 16[block 105]. In one embodiment, as indicated in [block 105] of FIG. 1A,point P may be at the approximate medial-lateral and anterior-posteriorcenter of the patient's joint 14. In other embodiments, point P may beat any other location in the 2D images 16, including anywhere on, nearor away from the bones 18, 20 or the joint 14 formed by the bones 18,20.

As described later in this overview, point P may be used to locate thecomputer generated 3D models 22, 28, 36 created from the 2D images 16and to integrate information generated via the 3D models. Depending onthe embodiment, point P, which serves as a position and/or orientationreference, may be a single point, two points, three points, a point plusa plane, a vector, etc., so long as the reference P can be used toposition and/or orient the 3D models 22, 28, 36 generated via the 2Dimages 16.

As discussed in detail below, the 2D images 16 are segmented along boneboundaries to create bone contour lines. Also, the 2D images 16 aresegmented along bone and cartilage boundaries to create bone andcartilage lines.

As shown in FIG. 1C, the segmented 2D images 16 (i.e., bone contourlines) are employed to create computer generated 3D bone-only (i.e.,“bone models”) 22 of the bones 18, 20 forming the patient's joint 14[block 110]. The bone models 22 are located such that point P is atcoordinates (X_(P), Y_(P), Z_(P)) relative to an origin (X₀, Y₀, Z₀) ofan X-Y-Z coordinate system [block 110]. The bone models 22 depict thebones 18, 20 in the present deteriorated condition with their respectivedegenerated joint surfaces 24, 26, which may be a result ofosteoarthritis, injury, a combination thereof, etc.

Computer programs for creating the 3D computer generated bone models 22from the segmented 2D images 16 (i.e., bone contour lines) include:Analyze from AnalyzeDirect, Inc., Overland Park, Kans.; Insight Toolkit,an open-source software available from the National Library of MedicineInsight Segmentation and Registration Toolkit (“ITK”), www.itk.org; 3DSlicer, an open-source software available from www.slicer.org; Mimicsfrom Materialise, Ann Arbor, Mich.; and Paraview available atwww.paraview.org. Further, some embodiments may use customized softwaresuch as OMSegmentation (renamed “PerForm” in later versions), developedby OtisMed, Inc. The OMSegmentation software may extensively use “ITK”and/or “VTK” (Visualization Toolkit from Kitware, Inc, available atwww.vtk.org.) Some embodiments may include using a prototype ofOMSegmentation, and as such may utilize InsightSNAP software.

As indicated in FIG. 1C, the 3D computer generated bone models 22 areutilized to create 3D computer generated “restored bone models” or“planning bone models” 28 wherein the degenerated surfaces 24, 26 aremodified or restored to approximately their respective conditions priorto degeneration [block 115]. Thus, the bones 18, 20 of the restored bonemodels 28 are reflected in approximately their condition prior todegeneration. The restored bone models 28 are located such that point Pis at coordinates (X_(P), Y_(P), Z_(P)) relative to the origin (X₀, Y₀,Z₀). Thus, the restored bone models 28 share the same orientation andpositioning relative to the origin (X₀, Y₀, Z₀) as the bone models 22.

In one embodiment, the restored bone models 28 are manually created fromthe bone models 22 by a person sitting in front of a computer 6 andvisually observing the bone models 22 and their degenerated surfaces 24,26 as 3D computer models on a computer screen 9. The person visuallyobserves the degenerated surfaces 24, 26 to determine how and to whatextent the degenerated surfaces 24, 26 surfaces on the 3D computer bonemodels 22 need to be modified to restore them to their pre-degeneratedcondition. By interacting with the computer controls 11, the person thenmanually manipulates the 3D degenerated surfaces 24, 26 via the 3Dmodeling computer program to restore the surfaces 24, 26 to a state theperson believes to represent the pre-degenerated condition. The resultof this manual restoration process is the computer generated 3D restoredbone models 28, wherein the surfaces 24′, 26′ are indicated in anon-degenerated state.

In one embodiment, the bone restoration process is generally orcompletely automated. In other words, a computer program may analyze thebone models 22 and their degenerated surfaces 24, 26 to determine howand to what extent the degenerated surfaces 24, 26 surfaces on the 3Dcomputer bone models 22 need to be modified to restore them to theirpre-degenerated condition. The computer program then manipulates the 3Ddegenerated surfaces 24, 26 to restore the surfaces 24, 26 to a stateintended to represent the pre-degenerated condition. The result of thisautomated restoration process is the computer generated 3D restored bonemodels 28, wherein the surfaces 24′, 26′ are indicated in anon-degenerated state. For more detail regarding a generally orcompletely automated system for the bone restoration process, see U.S.patent application Ser. No. 12/111,924 to Park, which is titled“Generation of a Computerized Bone Model Representative of aPre-Degenerated State and Usable in the Design and Manufacture ofArthroplasty Devices”, was filed Apr. 29, 2008, and is incorporated byreference in its entirety into this Detailed Description.

As depicted in FIG. 1C, the restored bone models 28 are employed in apre-operative planning (“POP”) procedure to determine saw cut locations30 and drill hole locations 32 in the patient's bones that will allowthe arthroplasty joint implants to generally restore the patient's jointline to it pre-degenerative alignment [block 120].

In one embodiment, the POP procedure is a manual process, whereincomputer generated 3D implant models 34 (e.g., femur and tibia implantsin the context of the joint being a knee) and restored bone models 28are manually manipulated relative to each other by a person sitting infront of a computer 6 and visually observing the implant models 34 andrestored bone models 28 on the computer screen 9 and manipulating themodels 28, 34 via the computer controls 11. By superimposing the implantmodels 34 over the restored bone models 28, or vice versa, the jointsurfaces of the implant models 34 can be aligned or caused to correspondwith the joint surfaces of the restored bone models 28. By causing thejoint surfaces of the models 28, 34 to so align, the implant models 34are positioned relative to the restored bone models 28 such that the sawcut locations 30 and drill hole locations 32 can be determined relativeto the restored bone models 28.

In one embodiment, the POP process is generally or completely automated.For example, a computer program may manipulate computer generated 3Dimplant models 34 (e.g., femur and tibia implants in the context of thejoint being a knee) and restored bone models or planning bone models 28relative to each other to determine the saw cut and drill hole locations30, 32 relative to the restored bone models 28. The implant models 34may be superimposed over the restored bone models 28, or vice versa. Inone embodiment, the implant models 34 are located at point P′ (X_(P′),Y_(P′), Z_(P′)) relative to the origin (X₀, Y₀, Z₀), and the restoredbone models 28 are located at point P (X_(P), Y_(P), Z_(P)). To causethe joint surfaces of the models 28, 34 to correspond, the computerprogram may move the restored bone models 28 from point P (X_(P), Y_(P),Z_(P)) to point P′ (X_(P′), Y_(P′), Z_(P′)), or vice versa. Once thejoint surfaces of the models 28, 34 are in close proximity, the jointsurfaces of the implant models 34 may be shape-matched to align orcorrespond with the joint surfaces of the restored bone models 28. Bycausing the joint surfaces of the models 28, 34 to so align, the implantmodels 34 are positioned relative to the restored bone models 28 suchthat the saw cut locations 30 and drill hole locations 32 can bedetermined relative to the restored bone models 28. For more detailregarding a generally or completely automated system for the POPprocess, see U.S. patent application Ser. No. 12/563,809 to Park, whichis titled Arthroplasty System and Related Methods, was filed Sep. 21,2009, and is incorporated by reference in its entirety into thisDetailed Description.

While the preceding discussion regarding the POP process is given in thecontext of the POP process employing the restored bone models ascomputer generated 3D bone models, in other embodiments, the POP processmay take place without having to employ the 3D restored bone models, butinstead taking placing as a series of 2D operations. For more detailregarding a generally or completely automated system for the POP processwherein the POP process does not employ the 3D restored bone models, butinstead utilizes a series of 2D operations, see U.S. patent applicationSer. No. 12/546,545 to Park, which is titled Arthroplasty System andRelated Methods, was filed Aug. 24, 2009, and is incorporated byreference in its entirety into this Detailed Description.

As indicated in FIG. 1E, in one embodiment, the data 44 regarding thesaw cut and drill hole locations 30, 32 relative to point P′ (X_(P′),Y_(P′), Z_(P′)) is packaged or consolidated as the “saw cut and drillhole data” 44 [block 145]. The “saw cut and drill hole data” 44 is thenused as discussed below with respect to [block 150] in FIG. 1E.

As can be understood from FIG. 1D, the 2D images 16 employed to generatethe bone models 22 discussed above with respect to [block 110] of FIG.1C are also segmented along bone and cartilage boundaries to form boneand cartilage contour lines that are used to create computer generated3D bone and cartilage models (i.e., “arthritic models”) 36 of the bones18, 20 forming the patient's joint 14 [block 130]. Like theabove-discussed bone models 22, the arthritic models 36 are located suchthat point P is at coordinates (X_(P), Y_(P), Z_(P)) relative to theorigin (X₀, Y₀, Z₀) of the X-Y-Z axis [block 130]. Thus, the bone andarthritic models 22, 36 share the same location and orientation relativeto the origin (X₀, Y₀, Z₀). This position/orientation relationship isgenerally maintained throughout the process discussed with respect toFIGS. 1B-1E. Accordingly, movements relative to the origin (X₀, Y₀, Z₀)of the bone models 22 and the various descendants thereof (i.e., therestored bone models 28, bone cut locations 30 and drill hole locations32) are also applied to the arthritic models 36 and the variousdescendants thereof (i.e., the jig models 38). Maintaining theposition/orientation relationship between the bone models 22 andarthritic models 36 and their respective descendants allows the “saw cutand drill hole data” 44 to be integrated into the “jig data” 46 to formthe “integrated jig data” 48 employed by the CNC machine 10 tomanufacture the customized arthroplasty jigs 2.

Computer programs for creating the 3D computer generated arthriticmodels 36 from the segmented 2D images 16 (i.e., bone and cartilagecontour lines) include: Analyze from AnalyzeDirect, Inc., Overland Park,Kans.; Insight Toolkit, an open-source software available from theNational Library of Medicine Insight Segmentation and RegistrationToolkit (“ITK”), www.itk.org; 3D Slicer, an open-source softwareavailable from www.slicer.org; Mimics from Materialise, Ann Arbor,Mich.; and Paraview available at www.paraview.org. Some embodiments mayuse customized software such as OMSegmentation (renamed “PerForm” inlater versions), developed by OtisMed, Inc. The OMSegmentation softwaremay extensively use “ITK” and/or “VTK” (Visualization Toolkit fromKitware, Inc, available at www.vtk.org.). Also, some embodiments mayinclude using a prototype of OMSegmentation, and as such may utilizeInsightSNAP software.

Similar to the bone models 22, the arthritic models 36 depict the bones18, 20 in the present deteriorated condition with their respectivedegenerated joint surfaces 24, 26, which may be a result ofosteoarthritis, injury, a combination thereof, etc. However, unlike thebone models 22, the arthritic models 36 are not bone-only models, butinclude cartilage in addition to bone. Accordingly, the arthritic models36 depict the arthroplasty target areas 42 generally as they will existwhen the customized arthroplasty jigs 2 matingly receive thearthroplasty target areas 42 during the arthroplasty surgical procedure.

As indicated in FIG. 1D and already mentioned above, to coordinate thepositions/orientations of the bone and arthritic models 36, 36 and theirrespective descendants, any movement of the restored bone models 28 frompoint P to point P′ is tracked to cause a generally identicaldisplacement for the “arthritic models” 36 [block 135].

As depicted in FIG. 1D, computer generated 3D surface models 40 of thearthroplasty target areas 42 of the arthritic models 36 are importedinto computer generated 3D arthroplasty jig models 38 [block 140]. Thus,the jig models 38 are configured or indexed to matingly receive thearthroplasty target areas 42 of the arthritic models 36. Jigs 2manufactured to match such jig models 38 will then matingly receive thearthroplasty target areas of the actual joint bones during thearthroplasty surgical procedure.

In one embodiment, the procedure for indexing the jig models 38 to thearthroplasty target areas 42 is a manual process. The 3D computergenerated models 36, 38 are manually manipulated relative to each otherby a person sitting in front of a computer 6 and visually observing thejig models 38 and arthritic models 36 on the computer screen 9 andmanipulating the models 36, 38 by interacting with the computer controls11. In one embodiment, by superimposing the jig models 38 (e.g., femurand tibia arthroplasty jigs in the context of the joint being a knee)over the arthroplasty target areas 42 of the arthritic models 36, orvice versa, the surface models 40 of the arthroplasty target areas 42can be imported into the jig models 38, resulting in jig models 38indexed to matingly receive the arthroplasty target areas 42 of thearthritic models 36. Point P′ (X_(P′), Y_(P′), Z_(P′)) can also beimported into the jig models 38, resulting in jig models 38 positionedand oriented relative to point P′ (X_(P′), Y_(P′), Z_(P′)) to allowtheir integration with the bone cut and drill hole data 44 of [block125].

In one embodiment, the procedure for indexing the jig models 38 to thearthroplasty target areas 42 is generally or completely automated, asdisclosed in U.S. patent application Ser. No. 11/959,344 to Park, whichis entitled System and Method for Manufacturing Arthroplasty Jigs, wasfiled Dec. 18, 2007 and is incorporated by reference in its entiretyinto this Detailed Description. For example, a computer program maycreate 3D computer generated surface models 40 of the arthroplastytarget areas 42 of the arthritic models 36. The computer program maythen import the surface models 40 and point P′ (X_(P′), Y_(P′), Z_(P′))into the jig models 38, resulting in the jig models 38 being indexed tomatingly receive the arthroplasty target areas 42 of the arthriticmodels 36. The resulting jig models 38 are also positioned and orientedrelative to point P′ (X_(P′), Y_(P′), Z_(P′)) to allow their integrationwith the bone cut and drill hole data 44 of [block 125].

In one embodiment, the arthritic models 36 may be 3D volumetric modelsas generated from a closed-loop process. In other embodiments, thearthritic models 36 may be 3D surface models as generated from anopen-loop process.

As indicated in FIG. 1E, in one embodiment, the data regarding the jigmodels 38 and surface models 40 relative to point P′ (X_(P′), Y_(P′),Z_(P′)) is packaged or consolidated as the “jig data” 46 [block 145].The “jig data” 46 is then used as discussed below with respect to [block150] in FIG. 1E.

As can be understood from FIG. 1E, the “saw cut and drill hole data” 44is integrated with the “jig data” 46 to result in the “integrated jigdata” 48 [block 150]. As explained above, since the “saw cut and drillhole data” 44, “jig data” 46 and their various ancestors (e.g., models22, 28, 36, 38) are matched to each other for position and orientationrelative to point P and P′, the “saw cut and drill hole data” 44 isproperly positioned and oriented relative to the “jig data” 46 forproper integration into the “jig data” 46. The resulting “integrated jigdata” 48, when provided to the CNC machine 10, results in jigs 2: (1)configured to matingly receive the arthroplasty target areas of thepatient's bones; and (2) having cut slots and drill holes thatfacilitate preparing the arthroplasty target areas in a manner thatallows the arthroplasty joint implants to generally restore thepatient's joint line to its pre-degenerated state.

As can be understood from FIGS. 1A and 1E, the “integrated jig data” 44is transferred from the computer 6 to the CNC machine 10 [block 155].Jig blanks 50 are provided to the CNC machine 10 [block 160], and theCNC machine 10 employs the “integrated jig data” to machine thearthroplasty jigs 2 from the jig blanks 50.

For a discussion of example customized arthroplasty cutting jigs 2capable of being manufactured via the above-discussed process, referenceis made to FIGS. 1F-1I. While, as pointed out above, the above-discussedprocess may be employed to manufacture jigs 2 configured forarthroplasty procedures involving knees, elbows, ankles, wrists, hips,shoulders, vertebra interfaces, etc., the jig examples depicted in FIGS.1F-1I are for total knee replacement (“TKR”) or partial knee replacement(“PKR”) procedures. Thus, FIGS. 1F and 1G are, respectively, bottom andtop perspective views of an example customized arthroplasty femur jig2A, and FIGS. 1H and 1I are, respectively, bottom and top perspectiveviews of an example customized arthroplasty tibia jig 2B.

As indicated in FIGS. 1F and 1G, a femur arthroplasty jig 2A may includean interior side or portion 100 and an exterior side or portion 102.When the femur cutting jig 2A is used in a TKR or PKR procedure, theinterior side or portion 100 faces and matingly receives thearthroplasty target area 42 of the femur lower end, and the exteriorside or portion 102 is on the opposite side of the femur cutting jig 2Afrom the interior portion 100.

The interior portion 100 of the femur jig 2A is configured to match thesurface features of the damaged lower end (i.e., the arthroplasty targetarea 42) of the patient's femur 18. Thus, when the target area 42 isreceived in the interior portion 100 of the femur jig 2A during the TKRor PKR surgery, the surfaces of the target area 42 and the interiorportion 100 match.

The surface of the interior portion 100 of the femur cutting jig 2A ismachined or otherwise formed into a selected femur jig blank 50A and isbased or defined off of a 3D surface model 40 of a target area 42 of thedamaged lower end or target area 42 of the patient's femur 18.

As indicated in FIGS. 1H and 1I, a tibia arthroplasty jig 2B may includean interior side or portion 104 and an exterior side or portion 106.When the tibia cutting jig 2B is used in a TKR or PKR procedure, theinterior side or portion 104 faces and matingly receives thearthroplasty target area 42 of the tibia upper end, and the exteriorside or portion 106 is on the opposite side of the tibia cutting jig 2Bfrom the interior portion 104.

The interior portion 104 of the tibia jig 2B is configured to match thesurface features of the damaged upper end (i.e., the arthroplasty targetarea 42) of the patient's tibia 20. Thus, when the target area 42 isreceived in the interior portion 104 of the tibia jig 2B during the TKRor PKR surgery, the surfaces of the target area 42 and the interiorportion 104 match.

The surface of the interior portion 104 of the tibia cutting jig 2B ismachined or otherwise formed into a selected tibia jig blank 50B and isbased or defined off of a 3D surface model 40 of a target area 42 of thedamaged upper end or target area 42 of the patient's tibia 20.

b. Automatic Segmentation of Scanner Modality Image Data to Generate 3DSurface Model of a Patient's Bone

In one embodiment as mentioned above, the 2D images 16 of the patient'sjoint 14 as generated via the imaging system 8 (see FIG. 1A and [block100] of FIG. 1B) are segmented or, in other words, analyzed to identifythe contour lines of the bones and/or cartilage surfaces that are ofsignificance with respect to generating 3D models 22, 36, as discussedabove with respect to [blocks 110 and 130] of FIGS. 1C and 1D.Specifically, a variety of image segmentation processes may occur withrespect to the 2D images 16 and the data associated with such 2D images16 to identify contour lines that are then compiled into 3D bone models,such as bone models 22 and arthritic models 36. A variety of processesand methods for performing image segmentation are disclosed in theremainder of this Detailed Description.

The imager 8 typically generates a plurality of image slices 16 viarepetitive imaging operations. Depending on whether the imager 8 is aMRI or CT imager, each image slice will be a MRI or CT slice. As shownin FIG. 2A, the image slice may depict the cancellous bone 200, thecortical bone 202 surrounding the cancellous bone, and the articularcartilage lining portions of the cortical bone 202 of an object ofinterest of a joint, e.g., a femur 204 in a patient's knee joint 14. Theimage may further depict the cancellous bone 206, the cortical bone 208of another object of interest in the joint, e.g., a tibia 210 of theknee joint 14. In one embodiment, each image slice 16 may be atwo-millimeter 2D image slice.

One embodiment may automatically segment one or more features ofinterest (e.g., bones) present in MRI or CT scans of a patient joint,e.g., knee, hip, elbow, etc. A typical scan of a knee joint mayrepresent approximately a 100-millimeter by 150-millimeter by150-millimeter volume of the joint and may include about 40 to 80 slicestaken in sagittal planes. A sagittal plane is an imaginary plane thattravels from the top to the bottom of the object (e.g., the human body),dividing it into medial and lateral portions. It is to be appreciatedthat a large inter-slice spacing may result in voxels (volume elements)with aspect ratios of about one to seven between the resolution in thesagittal plane (e.g., the y z plane) and the resolution along the x axis(i.e., each scan slice lies in the yz plane with a fixed value of x).For example, a two-millimeter slice that is 150-millimeters by150-millimeters may be comprised of voxels that are approximately0.3-millimeter by 0.3-millimeter by 2-millimeters (for a 512 by 512image resolution in the sagittal plane).

In one embodiment, each slice may be a gray scale image with aresolution of 512 by 512 voxels where the voxel value represents thebrightness (intensity) of the voxel. The intensity may be stored as a16-bit integer resulting in an intensity range from 0 to 65,535, where 0may represent black and 65,535 may represent white. The intensity ofeach voxel typically represents the average intensity of the voxelvolume. Other embodiments may employ scans having higher or lowerresolutions in the sagittal plane, different inter-slice spacing, orimages where the intensity may be represented by a 24 bit vector (e.g.,eight bits each for a red component, green component and bluecomponent). Additionally, other embodiments may store intensity valuesas 32-bit signed integers or floating point values.

Typical MRI and CT scan data generally provide images where parts of abone boundary of interest may be well defined while other parts of thebone boundary may be difficult to determine due to voxel volumeaveraging, the presence of osteophyte growth, the presence of tissuehaving similar image intensities in neighboring areas to the object tobe segmented, amongst other things. Such poor definition of parts of thebone boundary in the images may cause traditional automated segmentationtechniques to fail. For example, FIG. 2A depicts regions 212 within aslice where an object boundary may not be visible due to neighboringtissue having about the same intensity as the feature of interest.Depicted in FIG. 2B are regions 214 that may be extended into the slicefrom adjacent slices due to a high voxel aspect ratio. Depicted in FIG.2C is a region 216 of the bone boundary 218 that may disappear or loseregularity when the bone boundary 218 is approximately tangent to theslice.

One embodiment may employ image segmentation techniques using a goldentemplate to segment bone boundaries and provide improved segmentationresults over traditional automated segmentation techniques. Suchtechniques may be used to segment an image when similarity betweenpixels within an object to be identified may not exist. That is, thepixels within a region to be segmented may not be similar with respectto some characteristic or computed property such as a color, intensityor texture that may be employed to associate similar pixels intoregions. Instead, a spatial relationship of the object with respect toother objects may be used to identify the object of interest. In oneembodiment, a 3D golden template of a feature of interest to besegmented may be used during the segmentation process to locate thetarget feature in a target scan. For example, when segmenting a scan ofa knee joint, a typical 3D image of a known good femur (referred to as agolden femur template) may be used to locate and outline (i.e., segment)a femur in a target scan.

Generally, much of the tissues surrounding the cancellous and corticalmatter of the bone to be segmented may vary from one MRI scan to anotherMRI scan. This may be due to disease and/or patient joint position(e.g., a patient may not be able to straighten the joint of interestbecause of pain). By using surrounding regions that have a stableconnection with the bone (e.g., the grown golden and boundary goldenregions of the template as described in more detail below), theregistration may be improved. Additionally, use of these regions allowsthe bone geometry of interest to be captured during the segmentationrather than other features not of interest. Further, the segmentationtakes advantage of the higher resolution of features of interest incertain directions of the scan data through the use of a combination of2D and 3D techniques, that selectively increases the precision of thesegmentation as described in more detail below with respect to refiningthe bone registration using an artificially generated image.

The segmentation method employed by one embodiment may accommodate avariety of intensity gradients across the scan data. FIGS. 3A-C depictintensity gradients (i.e., the intensity varies non-uniformly across theimage) in slices (an intensity gradient that is darker on the top andbottom as depicted in FIG. 3A, an intensity gradient that is darker onthe bottom as depicted in FIG. 3B, and an intensity gradient 220 that isbrighter on the sides as depicted in FIG. 3C) that may be segmented byone embodiment. Further, the embodiment generally does not requireapproximately constant noise in the slices to be segmented. Theembodiment may accommodate different noise levels, e.g., high noiselevels as depicted in FIG. 4A as well as low noise levels as depicted inFIG. 4B. The decreased sensitivity to intensity gradients and noiselevel typically is due to image registration techniques using a goldentemplate, allowing features of interest to be identified even though thefeature may include voxels with differing intensities and noise levels.

Segmentation generally refers to the process of partitioning a digitalimage into multiple regions (e.g., sets of pixels for a 2D image or setsof voxels in a 3D image). Segmentation may be used to locate features ofinterest (bones, cartilage, ligaments, etc.) and boundaries (lines,curves, etc. that represent the bone boundary or surface) in an image.In one embodiment, the output of the automatic segmentation of the scandata may be a set of images (scan slices 16) where each image 16includes a set of extracted closed contours representing bone outlinesthat identify respective bone location and shape for bones of interest(e.g., the shape and location of the tibia and femur in the scan data ofa knee joint). The generation of a 3D model correspondent to the aboveclosed contours may be additionally included into the segmentationprocess. The automatic or semi-automatic segmentation of a joint, usingimage slices 16 to create 3D models (e.g., bone models 22 and arthriticmodels 36) of the surface of the bones in the joint, may reduce the timerequired to manufacture customized arthroplasty cutting jigs 2. It is tobe appreciated that certain embodiments may generate open contours ofthe bone shapes of interest to further reduce time associated with theprocess.

In one embodiment, scan protocols may be chosen to provide gooddefinition in areas where precise geometry reconstruction is requiredand to provide lower definition in areas that are not as important forgeometry reconstruction. The automatic or semi-automatic imagesegmentation of one embodiment employs components whose parameters maybe tuned for the characteristics of the image modality used as input tothe automatic segmentation and for the features of the anatomicalstructure to be segmented, as described in more detail below.

In one embodiment, a General Electric 3T MRI scanner may be used toobtain the scan data. The scanner settings may be set as follows: pulsesequence: FRFSE-XL Sagittal PD; 3 Pane Locator—Scout Scan Thickness:4-millimeters; Imaging Options: TRF, Fast, FR; Gradient Mode: Whole; TE:approximately 31; TR: approximately 2100; Echo Train Length: 8;Bandwidth: 50 Hz; FOV: 16 centimeters, centered at the joint line; PhaseFOV: 0.8 or 0.9; Slice Thickness: 2 millimeters; Spacing: Interleave;Matrix: 384×192; NEX: 2; Frequency: SI; and Phase Correct: On. It is tobe appreciated that other scanners and settings may be used to generatethe scan data.

Typically, the voxel aspect ratio of the scan data is a function of howmany scan slices may be obtained while a patient remains immobile. Inone embodiment, a two-millimeter inter-slice spacing may be used duringa scan of a patient's knee joint. This inter-slice spacing providessufficient resolution for constructing 3D bone models of the patient'sknee joint, while allowing sufficiently rapid completion of scan beforethe patient moves.

FIG. 5 depicts a MRI scan slice that illustrates image regions wheregood definition may be needed during automatic segmentation of theimage. Typically, this may be areas where the bones come in contactduring knee motion, in the anterior shaft area next to the joint andareas located at about a 10- to 30-millimeter distance from the joint.Good definition may be needed in regions 230 of the tibia 232 andregions 234 of the femur 236. Regions 238 depict areas where the tibiais almost tangent to the slice and boundary information may be lost dueto voxel volume averaging.

Voxel volume averaging may occur during the data acquisition processwhen the voxel size is larger than a feature detail to be distinguished.For example, the detail may have a black intensity while the surroundingregion may have a white intensity. When the average of the contiguousdata enclosed in the voxel is taken, the average voxel intensity valuemay be gray. Thus, it may not be possible to determine in what part ofthe voxel the detail belongs.

Regions 240 depict areas where the interface between the cortical boneand cartilage is not clear (because the intensities are similar), orwhere the bone is damaged and may need to be restored, or regions wherethe interface between the cancellous bone and surrounding region may beunclear due to the presence of a disease formation (e.g., an osteophytegrowth which has an image intensity similar to the adjacent region).

FIG. 6 depicts a flowchart illustrating one method for automatic orsemi-automatic segmentation of Femur and Tibia Planning models of animage modality scan (e.g., an MRI scan) of a patient's knee joint.Initially, operation 250 obtains a scan of the patient's knee joint. Inone embodiment, the scan may include about 50 sagittal slices. Otherembodiments may use more or fewer slices. Each slice may be a gray scaleimage having a resolution of 512 by 512 voxels. The scan may representapproximately a 100-millimeter by 150-millimeter by 150-millimetervolume of the patient's knee. While the invention will be described foran MRI scan of a knee joint, this is by way of illustration and notlimitation. The invention may be used to segment other types of imagemodality scans such as computed tomography (CT) scans, ultrasound scans,positron emission tomography (PET) scans, etc., as well as other jointsincluding, but not limited to, hip joints, elbow joints, etc. Further,the resolution of each slice may be higher or lower and the images maybe in color rather than gray scale. It is to be appreciated thattransversal or coronal slices may be used in other embodiments.

After operation 250 obtains scan data (e.g., scan images 16) generatedby imager 8, operation 252 may be performed to segment the femur data ofthe scan data. During this operation, the femur may be located andspline curves 270 may be generated to outline the femur shape or contourlines in the scan slices, as depicted in FIGS. 7A-7K. It should beappreciated that one or more spline curves may be generated in eachslice to outline the femur contour depending on the shape and curvatureof the femur as well as the femur orientation relative to the slicedirection.

Next, in operation 254, a trained technician may verify that thecontours of the femur spline curves generated during operation 252follow the surface of the femur bone. The technician may determine thata spline curve does not follow the bone shape in a particular slice. Forexample, FIG. 8 depicts an automatically generated femur spline curve274. The technician may determine that the curve should be enlarged inthe lower left part 276 of the femur. There may be various reasons whythe technician may decide that the curve needs to be modified. Forexample, a technician may want to generate a pre-deteriorated boneshape, yet the bone may be worn out in this region and may needreconstruction. The technician may determine this by examining theoverall 3D shape of the segmented femur and also by comparing lateraland medial parts of the scan data. The segmented region of the slice maybe enlarged by dragging one or more control points 278 located on thespline curve 274 to adjust the curve to more closely follow the femurboundary as determined by the technician, as shown by adjusted curve280. The number of control points on a spline curve may be dependent onthe curve length and curvature variations. Typically, 10-25 controlpoints may be associated with a spline curve for spline modification.

Once the technician is satisfied with all of the femur spline curves inthe scan slices, operation 256 generates a watertight triangular meshgeometry from the femur segmentation that approximates the 3D surface ofthe femur. The mesh closely follows the femur spline curves 270 andsmoothly interpolates between them to generate a 3D surface model of thefemur. FIG. 9 depicts typical 3D mesh geometry 290 of a target femurgenerated by one embodiment. Such a 3D model may be a 3D surface modelor 3D volume model resulting from open-loop contour lines or closed loopcontour lines, respectively. In one embodiment, such a 3D model asdepicted in FIG. 9 may be a bone model 22 or an arthritic model 36.

After operation 256, operation 258 may be performed to segment the tibiadata in the scan data. During this operation, the tibia is located andspline curves may be generated to locate and outline the shape of thetibia found in the scan slices, as depicted by tibia spline curves 272in FIGS. 7A-7K. It should be appreciated that one or more spline curvesmay be generated in each slice to outline the tibia depending on theshape and curvature of the tibia as well as the tibia orientationrelative to the slice direction.

Next, in operation 260, the technician may verify the tibia splinecurves generated during operation 258. The technician may determine thata spline curve does not follow the tibia in a particular slice. Forexample, referring back to FIG. 8, an automatically generated tibiaspline curve 282 is depicted that may not follow the tibia in the rightpart of the tibia due to the presence of an osteophyte growth 284. Thepresence of the osteophyte growth 284 may be determined by examiningneighboring slices. In this case, the segmented region may be reduced bydragging one or more control points 286 located on the spline curve tomodify the tibia spline curve 282 to obtain the adjusted tibia splinecurve 288. As previously discussed, each spline curve may haveapproximately 10-25 control points depending on the length and curvaturevariation of the spline curve.

When the purpose of the segmentation is generating bone models that willbe shown to a surgeon in the images where they are overlapped byimplants, a technician will not need to restore the segmentation modelto its pre-deteriorated bone shape, and thus will not need to spend timeon adjusting splines to follow the pre-deteriorated bone shape. Alsothere is no need to get highly precise segmentation in the bone areasthat are to be replaced with implant. So there is no need to spend timeon adjusting the non-perfect curves in the “to be replaced” areas.

Once the technician is satisfied with all of the tibia spline curves inthe scan slices, operation 262 generates a watertight triangular meshgeometry from the tibia segmentation. The mesh closely follows thespline curves and smoothly interpolates between them to generate a 3Dsurface model of the tibia. FIG. 10 depicts a typical 3D mesh geometry292 of a target tibia generated by one embodiment. Such a 3D model maybe a 3D surface model or 3D volume model resulting from open-loopcontour lines or closed loop contour lines, respectively. In oneembodiment, such a 3D model as depicted in FIG. 10 may be a bone model22 or an arthritic model 36.

Because the objects to be located in the scan data typically cannot besegmented by grouping similar voxels into regions, a golden templaterepresentative of a typical size and shape of the feature of interestmay be employed during the segmentation process to locate the targetfeature of interest.

FIG. 11 depicts a flowchart illustrating one method for generating agolden template. The method will be described for generating a goldentemplate of a tibia by way of illustration and not limitation. Themethod may be used to generate golden templates of other bonesincluding, but not limited to a femur bone, a hip bone, etc.

Initially, operation 300 obtains a scan of a tibia that is not damagedor diseased. The appropriate tibia scan may be chosen by screeningmultiple MRI tibia scans to locate a MRI tibia scan having a tibia thatdoes not have damaged cancellous and cortical matter (i.e., no damage intibia regions that will be used as fixed images to locate acorresponding target tibia in a target scan during segmentation), whichhas good MRI image quality, and which has a relatively average shape,e.g., the shaft width relative to the largest part is not out ofproportion (which may be estimated by eye-balling the images). Thistibia scan data, referred to herein as a golden tibia scan, may be usedto create a golden tibia template. It is to be appreciated that severalMRI scans of a tibia (or other bone of interest) may be selected, atemplate generated for each scan, statistics gathered on the successrate when using each template to segment target MRI scans, and the onewith the highest success rate selected as the golden tibia template.

In other embodiments, a catalog of golden models may be generated forany given feature, with distinct variants of the feature depending onvarious patient attributes, such as (but not limited to) weight, height,race, gender, age, and diagnosed disease condition. The appropriategolden mesh would then be selected for each feature based on a givenpatient's characteristics.

Then, in operation 302 the tibia is segmented in each scan slice. Eachsegmentation region includes the cancellous matter 322 and corticalmatter 324 of the tibia, but excludes any cartilage matter to form agolden tibia region, outlined by a contour curve 320, as depicted inFIG. 12A.

Next, operation 304 generates a golden tibia mesh 340 from theaccumulated golden tibia contours of the image slices, as illustrated inFIG. 13A.

Next, operation 306 increases the segmented region in each slice bygrowing the region to include boundaries between the tibia and adjacentstructures where the contact area is generally relatively stable fromone MRI scan to another MRI scan. This grown region may be referred toherein as a grown golden tibia region, outlined by contour curve 328, asdepicted in FIG. 12A.

The grown golden region may be used to find the surface that separatesthe hard bone (cancellous and cortical) from the outside matter(cartilage, tendons, water, etc.). The changes in voxel intensities whengoing from inside the surface to outside of the surface may be used todefine the surface. The grown golden region may allow the registrationprocess to find intensity changes in the target scan that are similar tothe golden template intensity changes near the surface. Unfortunately,the golden segmentation region does not have stable intensity changes(e.g., near the articular surface) or may not have much of an intensitychange. Thus, the grown region typically does not include such regionsbecause they do not provide additional information and may slow down theregistration due to an increased number of points to be registered.

Finally, use of a grown golden region may increase the distance wherethe metric function detects a feature during the registration process.When local optimization is used, the registration may be moved in aparticular direction only when a small movement in that directionimproves the metric function. When a golden template feature is fartheraway from the corresponding target bone feature (e.g., when there is asignificant shape difference), the metric typically will not move towardthat feature. Use of the larger grown region may allow the metric todetect the feature and move toward it.

Next, operation 308 cuts off most of the inner part of the grown goldentibia region to obtain a boundary golden tibia region 330 depicted inFIG. 12A. The boundary golden tibia region 330 is bounded on the insideby contour curve 332 and the outside by contour curve 328.

The boundary region may be used to obtain a more precise registration ofthe target bone by using the interface from the cancellous bone to thecortical bone. This may be done so that intensity variations in otherareas (e.g., intensity variations deep inside the bone) that may movethe registration toward wrong features and decrease the precision oflocating the hard bone surface are not used during the registration.

Then, operation 310 applies Gaussian smoothing with a standard deviationof two pixels to every slice of the golden tibia scan. In oneembodiment, a vtkImageGaussianSmooth filter (part of VisualizationToolkit, a free open source software package) may be used to perform theGaussian smoothing by setting the parameter “Standard Deviation” to avalue of two.

Then, operation 312 generates an anchor segmentation. The anchorsegmentation typically follows the original segmentation where the tibiaboundary is well defined in most MRI scans. In areas where the tibiaboundary may be poorly defined, but where there is another well definedfeature close to the tibia boundary, the anchor segmentation may followthat feature instead. For example, in an area where a healthy bonenormally has cartilage, a damaged bone may or may not have cartilage. Ifcartilage is present in this damaged bone region, the bone boundaryseparates the dark cortical bone from the gray cartilage matter. Ifcartilage is not present in this area of the damaged bone, there may bewhite liquid matter next to the dark cortical bone or there may beanother dark cortical bone next to the damaged bone area. Thus, theinterface from the cortical bone to the outside matter in this region ofthe damaged bone typically varies from MRI scan to MRI scan. In suchareas, the interface between the cortical and the inner cancellous bonemay be used. These curves may be smoothly connected together in theremaining tibia areas to obtain the tibia anchor segmentation curve 358,depicted in FIG. 14A.

Then, operation 314 may determine three disjoint regions along theanchor segmentation boundary. Each of these regions is generally welldefined in most MRI scans. FIG. 14A depicts these three disjoint regionsfor a particular image slice. The first region 350, referred to hereinas the tibia InDark-OutLight region, depicts a region where the anchorsegmentation boundary separates the inside dark intensity corticalmatter voxels from the outside light intensity voxels. The second region352, referred to herein as the tibia InLight-OutDark region, depicts aregion where the boundary separates the inside light intensitycancellous matter voxels from the outside dark intensity cortical mattervoxels. Finally, region 354, referred to herein as the tibiaDark-in-Light region, depicts a region that has a very thin layer ofdark intensity cortical matter voxels along the boundary, but which haslight intensity matter voxels away from the boundary (i.e., on bothsides of the boundary). Generally, the other regions along the anchorsegmentation boundary vary from scan to scan or may not be clear in mostof the scans, as depicted by regions 356. Such regions may be anosteophyte growth with an arbitrary shape but which has about the sameintensity as the region next to it. Thus, such regions typically are notused as anchor regions in one embodiment of the invention.

Finally, operation 316 generates a mesh corresponding to the anchorsegmentation and also generates a mesh for each anchor region. FIG. 15Adepicts the anchor segmentation mesh 360, the InDark-OutLight anchorregion mesh 362, the InLight-OutDark anchor region mesh 364 and theDark-in-Light anchor region mesh 366 for the tibia. These 3D meshesmodel the surface of the golden tibia in the specified regions. It is tobe appreciated that the 3D meshes are distinct and generally are notcombined to create a composite mesh. These meshes may be used to createan artificial fixed image that is used during the registration processas described in more detail below.

A golden template of a femur may also be generated in a similar mannerusing the method depicted by FIG. 11. FIG. 12B depicts the golden femurregion, outlined by a contour curve 320A, the grown femur region,outlined by contour curve 328A, and the boundary golden femur region330A bounded on the inside by contour curve 332A and the outside bycontour curve 328A. FIG. 13B depicts the golden femur mesh 340A. FIG.14B depicts the femur anchor segmentation curve 358A, the femurInDark-OutLight region 350A and the femur InLight-OutDark region 352A.Finally, FIG. 15B depicts the anchor segmentation mesh 360A, theInDark-OutLight anchor region mesh 362A and the InLight-OutDark anchorregion mesh 364A for the femur.

FIG. 16 depicts a flowchart illustrating one method for performingautomatic segmentation (e.g., operation 252 or operation 258 of FIG. 6)of the scan data of a joint (e.g., a MRI scan of a knee joint) usinggolden template registration. The segmentation method may be used tosegment the femur (operation 252 of FIG. 6) and/or the tibia (operation258 of FIG. 6) in either the left or right knee. Different goldentemplate data may be used to segment the left tibia, right tibia, leftfemur or right femur. Additionally, other embodiments may segment otherjoints, including but not limited to, hip joints, elbow joints, by usingan appropriate golden template of the feature of interest to besegmented.

Initially, operation 370 maps the segmented 3D golden template andmarked regions (e.g., grown and boundary regions) to the target scandata using image registration techniques. This may be done to locate thecorresponding feature of interest in the target scan (e.g., a targetfemur or tibia). Registration transforms the template image coordinatesystem into the target coordinate system. This allows the template imageto be compared and/or integrated with the target image.

Next, operation 372 refines the registration near the feature (e.g., abone) boundary of interest. Anchor segmentation and anchor regions maybe used with a subset of 3D free-form deformations to move points withinthe plane of the slices (e.g., the yz plane) but not transversal (alongthe x axis) to the slices. Refinement of the initial registrationoperation may be necessary to correct errors caused by a high voxelaspect ratio. When a point from a golden template is mapped onto thetarget scan, it generally maps to a point between adjacent slices of thescan data. For example, if a translation occurs along the x direction,then the point being mapped may only align with a slice when thetranslation is a multiple of the inter-slice scan distance (e.g., amultiple of two-millimeters for an inter-slice spacing oftwo-millimeters). Otherwise, the point will be mapped to a point thatfalls between slices. In such cases, the intensity of the target scanpoint may be determined by averaging the intensities of correspondingpoints (voxels) in the two adjacent slices. This may further reduceimage resolution. Additionally, refinement of the initial registrationoperation may correct for errors due to unhealthy areas and/or limitedcontrast areas. That is, the golden template may be partially pulledaway from the actual bone boundary in diseased areas and/or minimalcontrast areas (e.g., toward a diseased area having a differentcontrast) during the initial registration operation.

Next, operation 374 generates a polygon mesh representation of thesegmented scan data. A polygon mesh typically is a collection ofvertices, edges, and faces that may define the surface of a 3D object.The faces may consist of triangles, quadrilaterals or other simpleconvex polygons. In one embodiment, a polygon mesh may be generated byapplying the registration transform found during operation 372 to allthe vertices of a triangle golden template mesh (i.e., the surface ofthe mesh is composed of triangular faces). It is to be appreciated thatthe cumulative registration transform typically represents the transformthat maps the golden template into the target MRI scan with minimalmisalignment error.

Finally, operation 376 generates spline curves that approximate theintersection of the mesh generated by operation 374 with the target MRIslices. Note that these spline curves may be verified by the technician(during operation 254 or operation 260 of FIG. 6).

FIG. 17 depicts a flowchart illustrating one method for mapping thesegmented golden femur template regions into the target scan using imageregistration techniques. Registration may be thought of as anoptimization problem with a goal of finding a spatial mapping thataligns a fixed image with a target image. Generally several registrationoperations may be performed, first starting with a coarse imageapproximation and a low-dimensional transformation group to find a roughapproximation of the actual femur location and shape. This may be doneto reduce the chance of finding wrong features instead of the femur ofinterest. For example, if a free-form deformation registration wasinitially used to register the golden femur template to the target scandata, the template might be registered to the wrong feature, e.g., to atibia rather than the femur of interest. A coarse registration may alsobe performed in less time than a fine registration, thereby reducing theoverall time required to perform the registration. Once the femur hasbeen approximately located using a coarse registration, finerregistration operations may be performed to more accurately determinethe femur location and shape. By using the femur approximationdetermined by the prior registration operation as the initialapproximation of the femur in the next registration operation, the nextregistration operation may find a solution in less time.

In one embodiment, each registration operation may employ a registrationframework 390 as depicted in FIG. 18. The registration framework 390 mayemploy an image similarity-based method. Such a method generallyincludes a transformation model T(X) 392, which may be applied tocoordinates of a fixed (or reference) image 394 (e.g., a golden femurtemplate) to locate their corresponding coordinates in a target image396 space (e.g., a MRI scan), an image similarity metric 398, whichquantifies the degree of correspondence between features in both imagespaces achieved by a given transformation, and an optimizer 400, whichtries to maximize image similarity (or minimize an opposite function) bychanging the parameters of the transformation model 392. An interpolator402 may be used to evaluate target image intensities at non-gridlocations (e.g., reference image points that are mapped to target imagepoints that lie between slices). Thus, a registration frameworktypically includes two input images, a transform, a metric, aninterpolator and an optimizer.

Referring again to FIG. 17, operation 380 may approximately register agrown femur region in a MRI scan using a coarse registrationtransformation. In one embodiment, this may be done by performing anexhaustive translation transform search on the MRI scan data to identifythe appropriate translation transform parameters that minimizestranslation misalignment of the reference image femur mapped onto thetarget femur of the target image. This coarse registration operationtypically determines an approximate femur position in the MRI scan.

A translational transform, translates (or shifts) all elements of animage by the same 3D vector. That is, the reference femur may be mappedinto the target image space by shifting the reference femur along one ormore axes in the target image space to minimize misalignment. Duringthis operation the reference femur is not rotated, scaled or deformed.In one embodiment, three parameters for the translation transformationmay be generated: one parameter for each dimension that specifies thetranslation for that dimension. The final parameters of the translationtransform minimizing the misalignment of the mapped reference femurimage coordinates into the target image space may be stored.

Next, operation 382 further refines the image registration determined byoperation 380. This may be done by approximately registering the grownfemur region of the reference golden template femur into the target MRIscan data using a similarity transformation. In one embodiment, asimilarity transformation may be performed in 3D space. The referencegolden femur region may be rotated in 3D, translated in 3D andhomogeneously scaled to map its coordinates into the target MRI scandata to minimize misalignment between the reference golden femur regionand the corresponding region in the target MRI scan. In someembodiments, a center of rotation may be specified so that both therotation and scaling operations are performed with respect to thespecified center of rotation. In one embodiment, a 3D similaritytransformation, specified by seven parameters, may be used. Oneparameter specifies the scaling factor, three parameters specify aversor that represents the 3D rotation and three parameters specify avector that represents the 3D translation in each dimension. A versor isa unit quaternion that provides a convenient mathematical notation forrepresenting orientations and rotations of objects in three dimensions.

In one embodiment, local minimization techniques may be employed withthe similarity transformation to obtain a refined registration of thereference golden femur region onto the target MRI scan that is not farfrom the registration of the reference golden femur region onto thetarget MRI scan found in the previous operation 190 and used as theinitial starting approximation. Registering the grown golden femurregion may increase the distance where the metric function detects afeature during the registration process. When local optimization isused, the registration may be moved in a particular direction only whena small movement in that direction improves the metric function. When agolden femur template feature is farther away from the correspondingtarget femur feature (e.g., when there is a significant shapedifference), the metric typically will not move toward that feature. Useof the larger grown femur region may allow the metric to detect thefeature and move toward it.

After operation 382, operation 384 further refines the imageregistration of the golden femur into the target scan. In oneembodiment, an affine transformation may be used to register coordinatesof a boundary golden femur region of a golden femur template into thetarget MRI scan data. In one embodiment, the approximate femurregistration found during operation 382 may be used as the initialstarting approximation for the affine transformation.

An affine transformation typically is a linear transformation followedby a translation. The affine transformation preserves collinearitybetween points (i.e., three points which lie on a line continue to becollinear after the transformation) and ratios of distances along aline. In one embodiment, a 3D affine transformation, specified by 12parameters, may be utilized. Nine parameters of the affinetransformation specify the linear transformation (which may berepresented by a three by three matrix) and three parameters of theaffine transformation specify the 3D translation in each dimension. Theparameters of the affine transform that minimizes the misalignment ofthe boundary golden femur region mapped into the target MRI scan datamay be stored.

Finally, operation 386 further refines the image registration of theboundary golden femur region. In one embodiment, a spline transformationmay be used to register the coordinates of the boundary golden femurregion into the MRI scan data (target image space). In one embodiment, a3D B-Spline deformable transformation may be employed and thetransformation found in operation 384 may be used as the initialtransformation values for the 3D B-Spline deformable transformation.

A B-Spline deformable transformation typically is a free formdeformation of an object using a deformation field where a deformationvector is assigned to every point in space. For example, a 3D B-splinedeformable transform T may specify a 3D vector V(P) for every point P inthe original 3D space that is moved by T such that T:P→P+V(P).

In one embodiment, a B-Spline transformation may be specified with M×Nparameters, where M is the number of nodes in the B-Spline grid and N isthe dimension of the space. In one embodiment, a 3D B-Spline deformabletransformation of order three may be used to map every reference image3D point into the target MRI scan by a different 3D vector. The field ofthe vectors may be modeled using B-splines. Typically a grid J×K×L ofcontrol points may be specified where J, K, and L are parameters of thetransformation.

In one embodiment, splines of order three may be used with a grid 9×6×6of control points. That is, the transformation employs nine controlpoints in the medial/lateral direction (i.e., the x direction), and sixcontrol points in the other directions (i.e., y and z directions). Twocontrol points in each dimension (i.e., 2 of 9 in the x direction, 2 of6 in the y direction and 2 of 6 in the z direction) may be used tospecify boundary conditions. As such, the inner spline nodes may form agrid of size 7 by 4 by 4 and the boundary conditions increase the gridto size 9 by 6 by 6. The parametric set for this transformation has adimension of 3×9×6×6=972 (i.e., each dimension may have a 9×6×6 grid ofcontrol points). The final parameters of the spline transformation thatminimizes the misalignment between the reference golden femur templateand the target MRI scan data may be stored. This may be referred to asthe cumulative femur registration transform herein.

FIG. 19 depicts a flowchart illustrating one method for mapping thesegmented golden tibia template regions into the target scan using imageregistration techniques. Generally several registration operations maybe performed, first starting with a coarse image approximation and alow-dimensional transformation group to find a rough approximation ofthe actual tibia location and shape. This may be done to reduce thechance of finding wrong features instead of the tibia of interest. Forexample, if a free-form deformation registration was initially used toregister the golden tibia template to the target scan data, the templatemight be registered to the wrong feature, e.g., to a femur rather thanthe tibia of interest. A coarse registration may also be performed inless time than a fine registration, thereby reducing the overall timerequired to perform the registration. Once the tibia has beenapproximately located using a coarse registration, finer registrationoperations may be performed to more accurately determine the tibialocation and shape. By using the tibia approximation determined by theprior registration operation as the initial approximation of the tibiain the next registration operation, the next registration operation mayfind a solution in less time.

In one embodiment, each registration operation may employ a registrationframework 390 as depicted in FIG. 18. The registration framework 390 mayemploy an image similarity-based method. Such a method generallyincludes a transformation model T(X) 392, which may be applied tocoordinates of a fixed (or reference) image 394 (e.g., a golden tibiatemplate) to locate their corresponding coordinates in a target image396 space (e.g., a MRI scan), an image similarity metric 398, whichquantifies the degree of correspondence between features in both imagespaces achieved by a given transformation, and an optimizer 400, whichtries to maximize image similarity by changing the parameters of thetransformation model 392. An interpolator 402 may be used to evaluatetarget image intensities at non-grid locations (i.e., reference imagepoints that are mapped to target image points that lie between slices).Thus, a registration framework typically includes two input images, atransform, a metric, an interpolator and an optimizer.

The automatic segmentation registration process will be described usingscan data that includes a right tibia bone. This is by way ofillustration and not limitation. Referring again to FIG. 19, operation410 may approximately register a grown tibia region in a MRI scan usinga coarse registration transformation. In one embodiment, this may bedone by performing an exhaustive translation transform search on the MRIscan data to identify the appropriate translation transform parametersthat minimizes translation misalignment of the reference image tibiamapped onto the target tibia of the target image. This coarseregistration operation typically determines an approximate tibiaposition in the MRI scan. During this operation, the tibia of thereference image may be overlapped with the target tibia of the targetimage using a translation transformation to minimize translationalmisalignment of the tibias.

A translational transform, translates (or shifts) an image by the same3D vector. That is, the reference tibia may be mapped into the targetimage space by shifting the reference tibia along one or more axes inthe target image space to minimize misalignment. During this operationthe reference tibia is not rotated, scaled or deformed. In oneembodiment, three parameters for the translation transformation may begenerated, one parameter for each dimension that specifies thetranslation for that dimension. The final parameters of the translationtransform minimizing the misalignment of the mapped reference tibiaimage coordinates into the target image space may be stored.

Next, operation 412 further refines the image registration determined byoperation 410. This may be done by approximately registering the growntibia region of the reference golden tibia template into the target MRIscan data using a similarity transformation. In one embodiment, asimilarity transformation may be performed in 3D space. The referencegolden tibia region may be rotated in 3D, translated in 3D andhomogeneously scaled to map its coordinates into the target MRI scandata to minimize misalignment between the reference golden tibia regionand the corresponding region in the target MRI scan. In someembodiments, a center of rotation may be specified so that both therotation and scaling operations are performed with respect to thespecified center of rotation. In one embodiment, a 3D similaritytransformation, specified by seven parameters, may be used. Oneparameter specifies the scaling factor, three parameters specify aversor that represents the 3D rotation and three parameters specify avector that represents the 3D translation in each dimension. A versor isa unit quaternion that provides a convenient mathematical notation forrepresenting orientations and rotations of objects in three dimensions.

In one embodiment, local minimization techniques may be employed withthe similarity transformation to obtain a refined registration of thereference golden tibia region onto the target MRI scan that is not farfrom the registration of the reference golden tibia region onto thetarget MRI scan found in the previous operation 410 and used as theinitial starting approximation. Registering the grown golden tibiaregion may increase the distance where the metric function detects afeature during the registration process. When local optimization isused, the registration may be moved in a particular direction only whena small movement in that direction improves the metric function. When agolden tibia template feature is farther away from the correspondingtarget tibia feature (e.g., when there is a significant shapedifference), the metric typically will not move toward that feature. Useof the larger grown tibia region may allow the metric to detect thefeature and move toward it.

After operation 412, operation 414 further refines the imageregistration. In one embodiment, an affine transformation may be used toregister coordinates of a boundary golden tibia region of a golden tibiatemplate into the target MRI scan data. In one embodiment, theapproximate tibia registration found during operation 412 may be used asthe initial starting approximation for the affine transformation.

An affine transformation typically is a linear transformation followedby a translation. The affine transformation preserves collinearitybetween points (i.e., three points which lie on a line continue to becollinear after the transformation) and ratios of distances along aline. In one embodiment, a 3D affine transformation, specified by 12parameters, may be utilized. Nine parameters of the affinetransformation specify the linear transformation (which may berepresented by a three by three matrix) and three parameters of theaffine transformation specify the 3D translation in each dimension. Theparameters of the affine transform that minimizes the misalignment ofthe boundary golden tibia region mapped into the target MRI scan datamay be stored.

Finally, operation 416 further refines the image registration of theboundary golden tibia region. In one embodiment, a spline transformationmay be used to register the coordinates of the boundary golden tibiaregion into the MRI scan data (target image space). In one embodiment, a3D B-Spline deformable transformation may be employed and thetransformation found in operation 414 may be used as the initialtransformation values for the 3D B-Spline deformable transformation.

A B-Spline deformable transformation typically is a free formdeformation of an object using a deformation field where a deformationvector is assigned to every point in space. In one embodiment, aB-Spline transformation may be specified with M×N parameters, where M isthe number of nodes in the B-Spline grid and N is the dimension of thespace. In one embodiment, a 3D B-Spline deformable transformation oforder three may be used to map every reference image 3D point into thetarget MRI scan by a different 3D vector. The field of the vectors maybe modeled using B-splines. Typically a grid J×K×L of control points maybe specified where J, K, and L are parameters of the transformation.

In one embodiment, splines of order three may be used with a grid 9×6×6of control points. That is, the transformation employs nine controlpoints in the medial/lateral direction (i.e., the x direction, and sixcontrol points in the other directions (i.e., the y and z directions).Two control points in each dimension (i.e., 2 of 9 in the x direction, 2of 6 in the y direction and 2 of 6 in the z direction) may be used tospecify boundary conditions. As such, the inner spline nodes may form agrid of size 7 by 4 by 4 and the boundary conditions increase the gridto size 9 by 6 by 6. The parametric set for this transformation has adimension of 3×9×6×6=972. The final parameters of the splinetransformation that minimizes the misalignment between the referencegolden tibia template and the target MRI scan data may be stored. Thismay be referred to as the cumulative tibia registration transformherein.

The shape of the tibia may vary more from patient to patient than doesthe shape of the femur. As a result, the affine transformation may notprovide a close enough registration of the golden tibia template to thetarget tibia in the target scan. This may cause the Splinetransformation to find a local optimum that may be far from the actualtibia in some areas. In one embodiment, an additional registrationoperation between the affine transform and spline transform operationsmay be performed to more closely align the golden tibia and the targettibia, allowing the spline transform to converge to the correct localoptimum rather than a nearby (but wrong) local optimum.

The class of transforms utilized generally should allow more flexibility(or degrees of freedom) than the Affine transform and less flexibilitythan the B-spline transforms. The number of degrees of freedom generallyis equal to the number of transform parameters. In one embodiment, aclass of transforms with more than 12 parameters and less than 3×9×6×6parameters may be used. For example, a B-spline transform with fewercontrol points (than used in the subsequent spline transform) may beused for the additional transform operation. Alternatively, thedeformations may be modeled using quadric rather than cubic functions.

In another embodiment, several golden tibia templates may be used thatrepresent typical tibia variations, e.g., golden tibia templates forvarum, valgum and normal tibia. In one embodiment, each of the goldentibia templates may be used during the translation, similarity andaffine transform registration operations to find the template thatprovides the best match (e.g., best correlation) in the affine transformregistration operation. This template may then be used in the remainingregistration operations.

Finally, in one embodiment, the tibia registration may be improved byperforming the tibia segmentation after the femur segmentation andadding a restriction on the tibia registration transformations such thatthe tibia may not penetrate the femur. In one embodiment, this may beimplemented by introducing a penalty for the penetration. In the targetMRI all the voxels that lie inside the femur splines may be marked. Themetric functions, described in more detail below, that are used in theregistration operations may be modified to include a penalty term. Thepenalty term may be computed by selecting a set of points on theboundary of the golden template segmentation, applying a transform tothe set of points (in a similar way as the transform is applied to thesample points used in the correlation computations), determining if atransformed sample point falls into any of the marked voxels, and addinga large value to the penalty term for each transformed sample point thatfalls into any of the marked voxels.

In each of the above registration operations, a metric may be used toquantify the degree of correspondence between features in both thereference image and target image achieved by a given transformation. Inone embodiment, the metric quantitatively measures how well thetransformed golden template image fits the target image (e.g., a targetMRI scan) and may compare the gray-scale intensity of the images using aset of sample points in the golden template region to be registered.

FIG. 20 depicts a flowchart illustrating one method for computing themetric used by the registration operations described above. For aparticular registration operation, the metric may be computed in thesame way, but the metric may have different parameters specified for theparticular registration operation. The metric may be referred to hereinas “local correlation in sample points.” Initially, operation 420selects a set of sample points in the golden template region to beregistered.

For the translation and similarity transformations, the sample pointsmay be selected as follows. Initially, a rectilinear grid of L×M×N thatcovers the whole bone in 3D space may be used. L, M, and N may vary fromone to 16. In one embodiment, an eight by eight grid in every imageslice may be used to select uniform sample points in the grown goldenregion of the golden template. For each grid cell, the first samplepoint is selected. If the sample point falls within the grown goldenregion, it is used. If the sample point falls outside the golden region,it is discarded.

For the affine and spline transformations, the sample points may bedetermined by randomly selecting one out of every 32 points in theboundary golden region of the MRI slice.

Next, operation 422 groups the selected points into buckets. In oneembodiment, buckets may be formed as follows. First, the 3D space may besubdivided into cells using a rectilinear grid. Sample points thatbelong to the same cell are placed in the same bucket. It should benoted that sample points may be grouped into buckets to compensate fornon-uniform intensities in the MRI scan.

For example, MRI scan data may be brighter in the middle of the imageand darker towards the edges of the image. This brightness gradienttypically is different for different scanners and may also depend onother parameters including elapsed time since the scanner was lastcalibrated. Additionally, high aspect ratio voxels typically result invoxel volume averaging. That is, cortical bone may appear very dark inareas where its surface is almost perpendicular to the slice andgenerally will not be averaged with nearby tissues. However, corticalbone may appear as light gray in the areas where its surface is almosttangent to the slice and generally may be averaged with a large amountof nearby tissues.

Next, operation 424 sub-samples the target MRI slice. Sub-sampling thetarget space generally has the effect of smoothing the metric function.This may remove tiny local minima such that the local minimizationalgorithm converges to a deeper minimum. In one embodiment, duringoperations 410 and 412 (of FIG. 19), each slice may be sub-sampled withan eight by eight grid. During operations 414 and 416 (of FIG. 19), eachslice may be sub-sampled with a four by four grid. That is, during thesub-sampling, one point from every grid cell may be selected (e.g., thefirst point) and the remaining points in the grid cells may bediscarded.

Next, operation 426 computes a correlation of the intensities of thepoints in each bucket and their corresponding points in the target MRIscan (after mapping). The correlation (NC) metric may be expressed as:

${{NC}\left( {A,B} \right)} = {\frac{\sum\limits_{i = 1}^{N}{A_{i}B_{i}}}{\sqrt{\left( {\sum\limits_{i = 1}^{N}A_{i}^{2}} \right)\left( {\sum\limits_{i = 1}^{N}B_{i}^{2}} \right)}}\frac{{N{\sum{A_{i}B_{i}}}} - {\sum{A_{i}{\sum B_{i}}}}}{\sqrt{{N{\sum A_{i}^{2}}} - \left( {\sum A_{i}} \right)^{2}}\sqrt{{N{\sum B_{i}^{2}}} - \left( {\sum B_{i}} \right)^{2}}}}$

where A_(i) is the intensity in the i^(th) voxel of image A, B_(i) isthe intensity in the corresponding i^(th) voxel of image B and N is thenumber of voxels considered, and the sum is taken from i equals one toN. It should be appreciated that the metric may be optimal when imagedifferences are minimized (or when the correlation of image similaritiesis maximized). The NC metric generally is insensitive to intensityshifts and to multiplicative factors between the two images and mayproduce a cost function with sharp peaks and well defined minima.

Finally, operation 428 averages the correlations computed in everybucket with weights proportional to the number of sample points in thebucket.

It is to be appreciated that the above process for computing the metricmay compensate for non-uniform intensities, for example, those describedabove with respect to FIGS. 3A-3C, in the MRI scan data.

During the registration process, an optimizer may be used to maximizeimage similarity between the reference image and target image byadjusting the parameters of a given transformation model to adjust thelocation of reference image coordinates in the target image. In oneembodiment, the optimizer for a registration operation may use thetransformed image (e.g., the transformed golden template) from theprevious registration operation as its initial approximation. Then,local optimization techniques may be used to search for a local optimumnear the initial starting approximation. This may be done so that anypotential matches farther away from the feature of interest (e.g., thefemur or tibia in a knee joint) reliably found in an earlier operationmay be eliminated.

For the translation transformation, an exhaustive search may beperformed using a grid 10×10×10 of size 5-millimeter translationvectors. A translation for every vector in the grid may be performed andthe translation providing a maximum local correlation in sample pointsmay be selected as the optimum translation.

For the similarity transformation, a regular step gradient descentoptimizer may be used by one embodiment. A regular step gradient descentoptimizer typically advances transformation parameters in the directionof the gradient and a bipartition scheme may be used to compute the stepsize. The gradient of a function typically points in the direction ofthe greatest rate of change and whose magnitude is equal to the greatestrate of change.

For example, the gradient for a three dimensional space may be given by:

${\nabla{f\left( {x,y,z} \right)}} = {\left( {\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z}} \right).}$

That is, the gradient vector may be composed of partial derivatives ofthe metric function over all the parameters defining the transform. Inone embodiment the metric function may be a composition of an outer andN inner functions. The outer function may compute a metric valueaccording to operations 426 and 428 given the vectors {A_(i)} and{B_(i)}. The N inner functions may map N sample points from the fixed(reference) image A_(i) into the target image B_(i) using the transformand evaluate intensities of the target image B_(i) in the mapped points.Each of the inner functions generally depends on the transformparameters as well as on the point in the “from” space to which thetransform is applied. When computing the partial derivatives, the chainrule for computing a derivative of the function composition may be used.

To find a local minimum, parameter steps may be taken in the directionof the negative of the metric gradient (or the approximate gradient)over the transform parameter space at the current point. This generallyoptimizes the metric which typically has a local minimum when featuresof the reference image mapped into corresponding features of the targetimage have minimal misalignment).

The initial center of rotation for the similarity transformation (e.g.,operation 382 of FIG. 17) may be specified as the center of a boundingbox (or minimum sized cuboid with sides parallel to the coordinateplanes) that encloses the feature (e.g., a bone) registered in thetranslation registration (e.g., operation 380 of FIG. 17). Scalingcoefficients of approximately 40-millimeters may be used for the scalingparameters when bringing them together with translation parameters. Itis to be appreciated that the gradient computation generally relies on acustomized metric function in the parameter space, due to the fact thatwith a similarity transformation, the transform parameters do not havethe same dimensionality. The translation parameters have a dimension ofmillimeters, while the parameters for rotational angles and scaling donot have a dimension of millimeters. In one embodiment, a metric p maybe defined asρ=SQRT(X ² +Y ² +Z ²+(40-millimeter*A1)²+ . . . )

where X is the translation along the x axis, Y is the translation alongthe y axis, Z is the translation along the z axis, A₁ is the firstrotation angle, etc. A scaling coefficient of approximately40-millimeters may be used because it is approximately half the size ofthe bone (in the anterior/posterior and medial/lateral directions) ofinterest and results in a point being moved approximately 40-millimeterswhen performing a rotation of one radian angle.

In one embodiment, a maximum move of 1.5-millimeters may be specifiedfor every point, a relaxation factor may be set to 0.98 and a maximum of300 iterations may be performed to determine the parameters of thesimilarity transformation that results in minimal misalignment betweenthe reference image and target MRI scan.

For the affine transformation, a regular step gradient optimizer may beused by one embodiment. Scaling coefficients of approximately40-millimeters may be used for the matrix coefficients variations whenbringing them together with translation parameters. A maximum1.0-millimeter move for every point may be set for each iteration, therelaxation factor may be set to 0.98 and a maximum of 300 iterations maybe performed to determine the parameters of the affine transformationthat results in minimal misalignment.

For the B-spline transformation, a modified regular step gradientdescent optimizer may be used by one embodiment when searching for thebest B-spline deformable transformation. An MRI image gradient may oftenfollow the bone surface in diseased areas (e.g., where the bone contactsurface is severely damaged and/or where osteophytes have grown). Such agradient may cause deformations of the golden template that wouldintroduce large distortions in the segmented bone shape.

In one embodiment, the MRI image gradient may be corrected for suchdeformations by computing a normal to golden boundary vector field whereevery vector points towards the closest point in the golden templateshape found during the affine transformation (e.g., operation 384 ofFIG. 17). This may be done using a distance map (also referred to as adistance transform). A distance map supplies each voxel of the imagewith the distance to the nearest obstacle voxel (e.g., a boundary voxelin a binary image). In one embodiment, the gradient of the signeddistance map of the golden tibia region may be mapped using the affinetransformation found in operation 384 of FIG. 17. In one embodiment, asigned Danielsson distance map image filter algorithm may be used. Then,the MRI image gradient may be projected onto the vector field to obtainthe corrected gradient field. This corrected gradient field is parallelto the normal to golden boundary field and typically defines a very thinsubset of the set of B-spline transformations that may be spanned duringthe optimization.

Additionally, rather than computing one gradient vector for thetransform space and taking a step along it, a separate gradient may becomputed for every spline node. In one embodiment, order three B-splines(with J×K×L control nodes) may be used and J×K×L gradients may becomputed, one for each control point. At every iteration, each of thespline nodes may be moved along its respective gradient. This may allowthe spline curve to be moved in low contrast areas at the same time itis moved in high contrast areas. A relaxation factor of 0.95 may be usedfor each spline node. A maximum move of one-millimeter may be set forevery point during an iteration and a maximum of 20 iterations may beperformed to find the parameters of the B-spline transformation thatprovides minimal misalignment of the golden tibia region mapped into thetarget MRI scan.

Once the position and shape of the feature of interest of the joint hasbeen determined using image registration (operation 370 of FIG. 16), theregistration results may be refined using anchor segmentation and anchorregions (operation 372 of FIG. 16). FIG. 21 depicts a flowchartillustrating one method for refining the registration results usinganchor segmentation and anchor regions. Typically, during thisoperation, one more registration may be done using an artificiallygenerated image for the fixed image 394 of the registration framework390. Use of an artificial image may improve the overall segmentation byregistering known good regions that typically do not change from scan toscan to correct for any errors due to diseased and/or low contrast areasthat otherwise may distort the registration.

Additionally, the artificial image may be used to increase surfacedetection precision of articular surfaces and shaft middle regions. Theimage slices typically have higher resolution in two dimensions (e.g.,0.3-millimeter in the y and z dimensions) and lower resolution in thethird dimension (e.g., 2-millimeters in the x dimension). The articularsurfaces and shaft middle regions typically are well defined in theimage slices due to these surfaces generally being perpendicular to theslices. The surface detection precision may be improved using acombination of 2D and 3D techniques that preserves the in-sliceprecision by only moving points within slices rather than betweenslices. Further, a 3D B-spline transform may be used such that theslices are not deformed independently of one another. Since each slicemay not contain enough information, deforming each slice independentlymay result in the registration finding the wrong features. Instead, theslices as a whole may be deformed such that the registration remainsnear the desired feature. While each slice may be deformed differently,the difference in deformation between slices generally is small suchthat the changes from one slice to the next are gradual.

In one embodiment, the artificial image may comprise a set of dark andlight sample points that may be used by the metric. All dark points inthe artificial image may have the same intensity value (e.g., 100) andall light points in the artificial image may have the same intensityvalue (e.g., 200). It should be appreciated that the correlations aregenerally insensitive to scaling and zero shift. Thus, any intensityvalues may be used as long as the dark intensity value is less than thelight intensity value.

Initially, operation 430 may apply the cumulative registration transform(computed by operation 370 of FIG. 16) to an anchor segmentation meshand its three associated anchor region meshes (e.g., InDark-OutLightmesh, InLight-OutDark mesh and Dark-in-Light mesh) to generate atransformed anchor segmentation mesh and associated transformed anchorregion meshes (transformed InDark-OutLight anchor mesh, transformedInLight-OutDark anchor mesh and transformed Dark-in-Light anchor mesh)that lie in a space identical to the target image space.

Then, operation 432 generates random sample points lying within a thinvolume surrounding the transformed anchor segmentation mesh surface. Inone embodiment, this may be a volume having an outer boundary defined bythe anchor segmentation mesh surface plus 1.5-millimeters and an innerboundary defined by the anchor segmentation mesh surface minus1.5-millimeters, which may be referred to herein as the 1.5-millimeterneighborhood. The random sample points may be generated such that theyare within the image slices of the target scan but not between theslices. For example, the image slices may be transversal to the x-axiswith a spacing of 2-millimeters (at x-axis locations 0.0, 2.0, 4.0, . .. ). When a sample point is selected, its x-coordinate may be one of0.0, 2.0, 4.0, etc. but may not be 1.7, 3.0, or some non-multiple of2.0.

In one embodiment, voxels may be marked in every image slice that belongto the 1.5-millimeter neighborhood as follows. First, the intersectionof the transformed anchor mesh with every image slice may be found. Itshould be appreciated that the intersection of the anchor mesh with animage slice may be a polyline(s). Then, in each image slice, thepolyline segments may be traversed and all pixels that intersect withthe mesh may be marked. Next, a Dilate filter may be applied to themarked pixels of each image slice using a radius of 1.5-millimeters. TheDilate filter typically enlarges the marked region by adding all thepoints that lie within a 1.5-millimeter distance from the originallymarked points.

After operation 432, operation 434 determines if a sample point liesinside the transformed InDark-OutLight mesh surface. If operation 434determines that the sample point lies inside the transformedInDark-OutLight mesh surface, then operation 442 is performed. Ifoperation 434 determines that the sample point does not lie inside thetransformed InDark-OutLight mesh surface, then operation 436 isperformed.

Operation 442 determines if the sample point lies inside the transformedanchor segmentation mesh surface. If operation 442 determines that thesample point lies inside the transformed anchor segmentation meshsurface, then operation 446 is performed. If operation 442 determinesthat the sample point does not lie inside the transformed anchorsegmentation mesh surface, then operation 448 is performed.

Operation 436 determines if the sample point lies inside the transformedInLight-OutDark mesh surface. If operation 436 determines that thesample point lies inside the transformed InLight-OutDark mesh surface,then operation 444 is performed. If operation 436 determines that thesample point does not lie inside the transformed InLight-OutDark meshsurface, then operation 438 is performed.

Operation 444 determines if the sample point lies inside the transformedanchor segmentation mesh surface. If operation 444 determines that thesample point lies inside the transformed anchor segmentation meshsurface, then operation 448 is performed. If operation 444 determinessample point does not lie within the transformed anchor segmentationmesh surface, then operation 446 is performed.

Operation 438 determines if the sample point lies inside the transformedDark-In-Light mesh surface. If operation 438 determines that the samplepoint lies inside the transformed Dark-In-Light mesh surface, thenoperation 440 is performed. If operation 438 determines that the samplepoint does not lie inside the transformed Dark-In-Light mesh surface,then operation 450 is performed.

Operation 440 determines if the sample point is within 0.75-millimeterof the surface of the transformed anchor segmentation mesh. If operation440 determines that the sample point is within 0.75-millimeter of thesurface of the transformed anchor segmentation mesh, then operation 446is performed. If operation 440 determines that the sample point is notwithin 0.75-millimeter of the surface of the anchor segmentation mesh,then operation 450 is performed.

Operation 446 adds the sample point to the artificial image as a darkpoint. Then, operation 450 is performed.

Operation 448 adds the sample point to the artificial image as a lightsample point. Then, operation 450 is performed.

Operation 450 determines if there are more randomly generated samplespoints to be added to the artificial image. If operation 450 determinesthat there are more randomly generated sample points to be added to theartificial image, then operation 434 is performed. If operation 450determines that there are no more randomly generated sample points to beadded to the artificial image, then operation 452 is performed.

FIG. 22 depicts a set of randomly generated light sample points 460 anddark sample points 462 over the target MRI 464. In one embodiment,approximately 8,000 sample points (light and dark) may be generated overthe entire artificial image.

Referring again to FIG. 21, if operation 450 determines that there areno more randomly generated sample points to be added to the artificialimage, operation 452 registers the set of dark and light points to thetarget MRI scan. This operation may perform a registration similar tothe registration operation 196 (depicted in FIG. 17). In thistransformation, a subset of B-spline deformable transformations may beperformed that move points along their respective slices, but nottransversal to their respective slices.

In a B-spline deformable transform, a translation vector for everycontrol point (e.g., in the set of J×K×L control points) may bespecified. To specify a transform that moves any point in 3D space alongthe y and z slice coordinates but not along the x coordinate, arestriction on the choice of translation vectors in the control pointsmay be introduced. In one embodiment, only translation vectors with thex coordinate set equal to zero may be used to move points in the planeof the slice (e.g., the y and z directions) but not transversal to theslice (e.g., the x direction).

The use of anchor region meshes which typically are well pronounced inmost image scans may reduce registration errors due to unhealthy areasand/or areas with minimal contrast differences between the feature to besegmented and surrounding image areas. For example, in the area where ahealthy bone normally has cartilage, a damaged bone may or may not havecartilage. If cartilage is present in this damaged bone region, the boneboundary separates the dark cortical bone from the gray cartilagematter. If cartilage is not present in this area of the damaged bone,there may be white liquid matter next to the dark cortical bone or theremay be another dark cortical bone next to the damage bone area. Thus,the interface from the cortical bone to the outside matter in thisregion of the damaged bone typically varies from MRI scan to MRI scan.In such areas, the interface between the cortical and the innercancellous bone may be used as an anchor region.

The use of a subset of B-Spline deformable transforms may reduce errorsdue to the 2-millimeter spacing between image slices.

FIG. 23 depicts a flowchart illustrating one method for generatingspline curves outlining the surface of a feature of interest in eachtarget MRI slice (e.g., operation 376 of FIG. 16). Initially, operation470 intersects the generated 3D mesh model of the feature surface with aslice of the target scan data. The intersection defines a polyline curveof the surface of the feature (e.g., bone) in each slice. Two or morepolyline curves may be generated in a slice when the bone is not verystraightly positioned with respect to the slice direction.

A polyline curve is a piecewise linear approximation to a curved featureshape. Generally, this curve should be easy to manipulate with a set ofcontrol points. The polyline curve may have many segments, making itmore difficult to manipulate the polyline curve (e.g., during operation254 or 260 of FIG. 6). One embodiment may generate one or more Kochaneksplines from the polyline curve. Each spline typically has a smallernumber of control points and typically fits the polyline curve withabout 0.2-millimeter deviation. Generally, a Kochanek spline may havemore control points along the high curvature regions of the polylinecurve and fewer control points along low curvature regions (i.e., wherethe curve tends to be flatter) of the polyline curve.

Once a polyline curve has been generated, operation 472 may compute apolyline parameterization, L_(i), as a function of the polyline'slength. FIG. 24 depicts a polyline curve 481 with n vertices, V₀, V₁, .. . V_(i−1), V_(i) . . . V_(n-1). Note that vertex V₀ follows vertexV_(n-1) to form a closed contour curve. The length of a segmentconnecting vertices Vi−1 and Vi may be denoted by ΔL, such that thelength parameterization, L_(i), of the polyline at vertex V_(i) may beexpressed as:L _(i) =ΔL ₀ +ΔL ₁ + . . . +ΔL _(i).

Next, operation 474 may compute a polyline parameterization, A_(i), as afunction of the polyline's tangent variation. The absolute value of theangle between a vector connecting vertices V_(i−1) and V_(i) and avector connecting vertices V_(i) and V_(i+1) may be denoted by ΔA_(i)such that the tangent variation parameter A_(i) at vertex V_(i) may beexpressed as:A _(i) =ΔA ₀ +ΔA ₁ + . . . +ΔA _(i).

Then, operation 476 determines a weighted sum parameterization of thepolyline length and tangent variation parameterizations. In oneembodiment the weighted sum parameterization, W_(i), at vertex V_(i) maybe computed as:W _(i) =α*L _(i) +β*A _(i)

where α may be set to 0.2 and β may be set to 0.8 in one embodiment.

Then, operation 478 may perform a uniform sampling of the polyline usingthe W parameterization results determined by operation 476. In oneembodiment, a spacing interval of approximately 3.7 of the W parametervalue may be used for positioning K new sample points. First, K may becomputed as follows:K=ROUND(W _(n)/3.7).

That is, the W parameter value, which is the last computed value W_(n),may be divided by 3.7 and the result rounded to the nearest integer toget the number of new sample points. Then, the spacing of the samplepoints, ΔW may be computed as:ΔW=W _(n) /K.

Finally, the K new sample points, which are uniformly spaced, may bepositioned at intervals ΔW of the parameter W. The resulting samplepoints may be used as control points for the Kochanek splines to convertthe polyline into a spline. A Kochanek spline generally has a tension, abias and a continuity parameter that may be used to change the behaviorof the tangents. That is, a closed Kochanek spline with K control pointstypically is interpolated with K curve segments. Each segment has astarting point, an ending point, a starting tangent and an endingtangent. Generally, the tension parameter changes the length of thetangent vectors, the bias parameter changes the direction of the tangentvectors and the continuity parameter changes the sharpness in changebetween tangents. In certain embodiments, the tension, bias andcontinuity parameters may be set to zero to generate a Catmull-Romspline.

In one embodiment, operation 478 may perform a linear interpolation ofW_(i) and W_(i+1) to locate a sample point that lies between W_(i) andW_(i+1). The interpolated value of W may be used to determine thecorresponding sample location in the segment connecting vertices V_(i)and V_(i+1).

In certain embodiments, operation 478 may divide the W parameter valueby six to obtain the new number of sample points K. That is,K=ROUND(W _(n)/6).

Then, a measure of closeness (i.e., how closely the spline follows thepolyline) may be computed as follows. First, the spline is sampled suchthat there are seven sample points in every arc of the spline (i.e., 7*Ksample points). Then, the sum of the squared distances of the samplepoints to the polyline may be computed. Next, the coordinates of the Kcontrol points are varied (i.e., two*K parameters). Then, a localoptimization algorithm is used to find the closest spline. If theclosest spline found during the optimization is not within a certainprecision (e.g., within approximately 0.4-millimeter of the polyline),then the number of control points, K, may be increased by one. The newnumber of control points may be uniformly distributed along the Wparameter, and another optimization performed to find the new closestspline. Generally one to two optimizations provide a spline that followsthe polyline with the desired degree of precision (e.g., withinapproximately 0.2-millimeter).

Finally, operation 480 determines if a spline curve(s) should begenerated for another image slice. If operation 480 determines that aspline curve should be generated for another slice, then operation 472is performed. If operation 480 determines that there are no more imageslices to be processed, the method terminates.

As discussed above, in one embodiment, the output of the segmentationmay be a triangular mesh (e.g., a 3D surface model) of the segmentedbone(s) of a joint (e.g., the femur and tibia of a knee joint). The meshgenerated generally represents a watertight surface that closely followsthe segmentation contour curves of the slices, smoothly interpolatesbetween the segmentation contour curves, and may have a low triangularcount.

In one embodiment, a triangular mesh may be generated as follows. Thesegmentation data may be represented in 3D using (x, y, z) coordinateswith the image slices transversal to the x direction. Thus, thesegmentation contours lie in yz planes with fixed values of x.Initially, an in-slice distance image may be computed for each segmentedslice. The value of each (y, z) pixel in an in-slice distance image isthe distance to the closest point in the contours when the point islocated inside one of the contours and is the inverse (i.e., negative)of the distance to the closest point in the contours when the point isoutside all of the contours.

Then, a marching cubes algorithm may be applied to the in-slice distanceimages to generate the mesh. The marching cubes algorithm is a computeralgorithm for extracting a polygonal mesh of an isosurface (i.e., thecontours) from a three-dimensional scalar field (or voxels). Thealgorithm typically proceeds through the voxels, taking eight neighborvoxels at a time (thus forming an imaginary cube) and determines thepolygon(s) needed to represent the part of the isosurface (i.e.,contour) that passes through the imaginary cube. The individual polygonsare then fused into the desired surface. The generated mesh generallypasses through the zero level of the signed distance function in eachslice such that the mesh lies close to the contours.

It is to be appreciated that the image resolution in the y and zdirections typically determines how well the zero level of the signeddistance function approximates the original contours and may alsodetermine the triangular count in the resulting mesh. In one embodiment,a voxel size of 1.5-millimeters in the y and z directions may be used.This typically yields deviations within 0.1-millimeter of the originalcontours and produces a smooth mesh.

In one embodiment, a smoothing operation may be performed in the xdirection (i.e., transversal to the image slices) to compensate forsurface waviness that may have been introduced when the automaticallygenerated contours were adjusted (e.g., during operation 260 of FIG. 6).Such waviness may occur in regions of an image slice where there isminimal contrast variation and the curve is positioned by thetechnician. Typically a smooth best guess mesh in uncertain areas may bedesired when generating a planning model that may be used to locate theposition of an implant. Alternatively, a smooth overestimation may bedesired in uncertain areas such as in an arthritic model used to createa jig.

In one embodiment, simple smoothing may be used and the amount ofsmoothing (i.e., how much a voxel value may be modified) may becontrolled by two user specified parameters, MaxUp and MaxDown. After anaverage is computed for a voxel, it is clamped using these values tolimit the amount of smoothing. The smoothing operation typically doesnot change the image much in areas where the image contrast is good. Forsmooth best guess averaging in uncertain areas, MaxUp and MaxDown mayeach be set to 1 millimeter. For smooth overestimation averaging inuncertain regions, MaxUp may be set to 2-millimeters and MaxDown may beset to 0-millimeter.

The operation of adjusting segments of the segmentation process will nowbe described with reference to FIG. 25, which depicts a flowchart forone method of adjusting segments (e.g., operation 254 or operation 260of the flowchart depicted in FIG. 6). In one embodiment, thesegmentation data may be manually adjusted by a trained techniciansitting in front of a computer 6 and visually observing theautomatically generated contour curves in the image slices on a computerscreen 9. By interacting with computer controls 11, the trainedtechnician may manually manipulate the contour curves. The trainedtechnician may visually observe all of the contours as a 3D surfacemodel to select an image slice for further examination.

Initially, in operation 482 a slice is selected for verification. In oneembodiment, the slice may be manually selected by a technician.

Next, operation 484 determines if the segmentation contour curve in theselected slice is good. If operation 484 determines that thesegmentation contour curve is good, then operation 494 is performed. Ifoperation 484 determines that the segmentation contour curve is notgood, then operation 486 is performed.

Operation 486 determines if the segmentation contour curve isapproximately correct. If operation 486 determines that the contourcurve is approximately correct, then operation 492 is performed.

In operation 492 incorrect points of the segmentation contour curve maybe repositioned. In one embodiment this may be performed manually by atrained technician. It is to be appreciated that it may be difficult forthe technician to determine where the correct contour curve should belocated in a particular slice. This may be due to missing or unclearbone boundaries and/or areas with little contrast to distinguish imagefeatures. In one embodiment, a compare function may be provided to allowthe technician to visually compare the contour curve in the currentslice with the contour curves in adjacent slices. FIG. 26 depicts animage showing the contour curve 510 (e.g., a spline curve) with controlpoints 512 of the contour curve 510 for the current image slice as wellthe contour curves 514, 516 of the previous and next image slices,respectively, superimposed on the current image slice.

It may be difficult to determine where the correct segmentation contourcurve should be located due to missing or unclear bone boundaries due tothe presence of unhealthy areas, areas with limited contrastdifferences, and/or voxel volume averaging. When visually comparingadjacent slices, the technician may visualize the data in 2D planes (xy,yz, and xz) and in 3D. In one embodiment, the technician may select anarea for examination by positioning a crosshair on a location in anywindow and clicking a mouse button to select that image point. Thecrosshair will be placed at the desired point and may be used toindicate the same location when the data is visualized in each window.

The technician may use the spline control points to manipulate the shapeof the curve. This may be done by using a mouse to click on a controlpoint and dragging it to a desired location. Additionally, thetechnician may add or delete spline curve control points. This may bedone by using a mouse to select two existing control points betweenwhich a control point will be inserted or deleted. Alternatively, thetechnician may use a mouse cursor to point to the location on the curvewhere a control point is to be inserted. In one embodiment, by pressingthe letter I on a keyboard and then positioning the cursor at thedesired location, clicking the left mouse button will insert the controlpoint. A control point may be deleted by pressing the letter D on thekeyboard and then positioning the cursor over the desired control pointto be deleted. The selected control point will change color. Theselected control point will be deleted when the left mouse button isclicked.

Referring again to FIG. 25, if operation 486 determines that the contourcurve is not approximately correct, operation 488 is performed to deletethe curve. Then, operation 490 is performed.

Operation 490 generates a new segmentation contour curve for the imageslice. In one embodiment, a technician may use a spline draw tool toinsert a new spline curve. With the spline draw tool, the technician mayclick on consecutive points in the current slice to indicate where thespline curve should be located and a spline curve is generated thatpasses through all of the indicated points. A right mouse click may beused to connect the first and last points of the new spline curve.Alternatively, the technician may use a paste command to copy the splinecurve(s) from the previous slice into the current slice. The splinecontrol points may then be manipulated to adjust the spline curves tofollow the feature in the current image slice.

In another embodiment, a paste similar command may be used by thetechnician to copy the spline curve from the previous slice into thecurrent slice. Rather then pasting a copy of the spline curve from theprevious slice, the spline curve may be automatically modified to passthrough similar image features present in both slices. This may be doneby registering a region around the spline curve in the previous slicethat is from about 0.7-millimeter outside of the curve to about5.0-millimeter within the curve. Initially, this region is registeredusing an affine transformation. Then, the result of the affine transformmay be used as a starting value for a B-Spline deformabletransformation. The metric used for the transform may be the localcorrelation in sample points metric described previously. Typically,more sample points may be taken closer to the curve and fewer samplepoints taken farther away from the curve. Next, the spline controlpoints may be modified by applying the final transformation found to thespline control points. Additionally, the trained technician may adjustfrom zero to a few control points in areas where the bone boundarychanges a lot from the slice due to the bone being tangent to the sliceor in areas of limited contrast (e.g., where there is an osteophytegrowth). Then, operation 492 is performed.

Operation 494 determines if there are additional slices to be verified.If operation 494 determines that there are additional slices to beverified, operation 482 is performed.

If operation 494 determines that there are no more slices to beverified, then operation 496 is performed. Operation 496 generates a 3Dsurface model of the segmented bone.

Then, operation 498 determines if the 3D surface model is good. In oneembodiment, a technician may manually determine if the 3D surface modelis good. The technician may use a spline 3D visualization tool thatgenerates a slice visualization showing the voxels inside all of thesplines in 3D, as illustrated by the 3D shape 520 depicted in FIG. 27.This spline 3D visualization tool typically may be generated in realtime to provide interactive updates to the technician as the splinecurves are manually edited. Alternatively, a mesh visualization may begenerated in response to a technician command. The mesh visualizationtypically generates a smooth mesh that passes close to all the splinecurves, e.g., mesh 290 depicted in FIG. 9.

If operation 498 determines that the 3D model is not good, thenoperation 500 is performed. Operation 500 selects a slice lying in anarea where the 3D shape is not good. In one embodiment, a technician maymanually select the slice. Then, operation 482 is performed.

If operation 498 determines that the 3D model is good, then the methodterminates.

The 3D surface models of the lower end of the femur and the upper end ofthe tibia of a patient's knee may be used to create arthroplasty jigsand/or implants. For example, the models may be used to create femur andtibia jigs that can be used with a patient's femur and tibia asdisclosed in the various U.S. Patent Applications incorporated byreference herein in this Detailed Description and filed by Park and Parket al. Automatic segmentation of image data to generate 3D bone modelsmay reduce the overall time required to perform a reconstructive surgeryto repair a dysfunctional joint and may also provide improved patientoutcomes.

c. Segmentation Using Landmarks of Scanner Modality Image Data toGenerate 3D Surface Model of a Patient's Bone

Now begins a discussion of an alternative embodiment of imagesegmentation. The alternative embodiment includes placing landmarks 777(in FIG. 35A-FIG. 35H) on image contours. The landmarks 777 are thenused to modify a golden bone model (e.g., golden femur or golden tibia),the resulting modified golden bone model being the output ofsegmentation.

Similar to the embodiment of image segmentation discussed above insection b. of this Detailed Discussion, in one version of thealternative embodiment of image segmentation, the 2D images 16 of thepatient's joint 14 are generated via the imaging system 8 (see FIG. 1Aand [block 100] of FIG. 1B). These images 16 are analyzed to identifythe contour lines of the bones and/or cartilage surfaces that are ofsignificance with respect to generating 3D models 22, 36, as discussedabove in section a. of this Detailed Discussion with respect to [blocks110 and 130] of FIGS. 1C and 1D. Specifically, a variety of imagesegmentation processes may occur with respect to the 2D images 16 andthe data associated with such 2D images 16 to identify contour linesthat are then compiled into 3D bone models, such as bone models 22,restored bone models 28, and arthritic models 36.

Algorithms and software are described in section b. of this DetailedDiscussion for automatic and semi-automatic image segmentation. Insection c. of the Detailed Description, alternative software tools andunderlying methods are described, such alternative tools and methodshelping a user to quickly generate bone models. Because the alternativesoftware requires some user input such as, for example, initial Landmarkpositions, final verification and, in some instances, adjustment, thisalternative segmentation process can be considered a semi-automaticsegmentation process.

In some cases the alternative embodiment described in section c. of thisDetailed Discussion may significantly reduce the user time spent onsegmentation. In particular, compared to manual segmentation (where theuser draws contour(s) by hand on each applicable slice for eachapplicable bone), the user time may be reduced by approximately fivetimes when a user segments a planning model intended for communicating apreoperative planning model to a surgeon. For that purpose, a user maygenerate 3D bone models with high precision in particular areas and lessprecision in other areas. In some implementations, a user may get highprecision (e.g., 0.5 mm) at well defined bone contours in MRI imagesoutside the implant regions and less precision (e.g., up to 2 mm) in theregions that will be replaced with implants by spending approximately3-4 minutes in the user interface (“UI”) setting landmarks for thealgorithm. If improved precision is desired, the user may position morelandmarks and thus spend more time in the UI.

In one embodiment, the software tool described in section c. of theDetailed Discussion is called “Segmentation using Landmarks”. This toolmay be implemented inside software application PerForm 1.0. A variety ofprocesses and methods for performing image segmentation using landmarksare disclosed in the remainder of this Detailed Description.

The imager 8 typically generates a plurality of image slices 16 viarepetitive imaging operations. Depending on whether the imager 8 is aMRI or CT imager, each image slice will be a MRI or CT slice. As shownin FIG. 2A, the image slice may depict the cancellous bone 200, thecortical bone 202 surrounding the cancellous bone, and the articularcartilage lining portions of the cortical bone 202 of an object ofinterest of a joint, e.g., a femur 204 in a patient's knee joint 14. Theimage may further depict the cancellous bone 206, the cortical bone 208of another object of interest in the joint, e.g., a tibia 210 of theknee joint 14. In one embodiment, each image slice 16 may be atwo-millimeter 2D image slice.

One embodiment may segment one or more features of interest (e.g.,bones) present in MRI or CT scans of a patient joint, e.g., knee, hip,elbow, etc. A typical scan of a knee joint may represent approximately a100-millimeter by 150-millimeter by 150-millimeter volume of the jointand may include about 40 to 80 slices taken in sagittal planes. Asagittal plane is an imaginary plane that travels from the top to thebottom of the object (e.g., the human body), dividing it into medial andlateral portions. It is to be appreciated that a large inter-slicespacing may result in voxels (volume elements) with aspect ratios ofabout one to seven between the resolution in the sagittal plane (e.g.,the y z plane) and the resolution along the x axis (i.e., each scanslice lies in the yz plane with a fixed value of x). For example, atwo-millimeter slice that is 150-millimeters by 150-millimeters may becomprised of voxels that are approximately 0.3-millimeter by0.3-millimeter by 2-millimeters (for a 512 by 512 image resolution inthe sagittal plane).

In one embodiment, each slice may be a gray scale image with aresolution of 512 by 512 voxels where the voxel value represents thebrightness (intensity) of the voxel. The intensity may be stored as a16-bit integer resulting in an intensity range from 0 to 65,535, where 0may represent black and 65,535 may represent white. The intensity ofeach voxel typically represents the average intensity of the voxelvolume. Other embodiments may employ scans having higher or lowerresolutions in the sagittal plane, different inter-slice spacing, orimages where the intensity may be represented by a 24 bit vector (e.g.,eight bits each for a red component, green component and bluecomponent). Additionally, other embodiments may store intensity valuesas 8-bit or 32-bit signed or unsigned integers or floating point values.

Typical MRI and CT scan data generally provide images where parts of abone boundary of interest may be well defined while other parts of thebone boundary may be difficult to determine due to voxel volumeaveraging, the presence of osteophyte growth, the presence of tissuehaving similar image intensities in neighboring areas to the object tobe segmented, amongst other things. Such poor definition of parts of thebone boundary in the images may cause fully automated segmentationtechniques to fail. For example, FIG. 2A depicts regions 212 within aslice where an object boundary may not be visible due to neighboringtissue having about the same intensity as the feature of interest.Depicted in FIG. 2B are regions 214 that may be extended into the slicefrom adjacent slices due to a high voxel aspect ratio. Depicted in FIG.2C is a region 216 of the bone boundary 218 that may disappear or loseregularity when the bone boundary 218 is approximately tangent to theslice.

In one embodiment, a user may provide some additional input to theauto-segmentation algorithm, and the algorithm could use the additionaluser input for more accurate and faster segmentation of features ofinterest. For example, the additional user input may be a set of pointson the boundary of the feature of interest. In the context of a kneeprocedure, the points might be on the Femur knee bone boundary or on theTibia knee bone boundary. These can be called landmark points or simplylandmarks 777.

In order for a user to provide particular landmark points, the softwaremay allow loading MRI or CT image data, viewing and scrolling over imageslices, specifying landmark points in the slices and editing them. Thesoftware may also allow visualization of the segmentation results (i.e.,segmentation curves drawn in the image slices). The software may alsogenerate a 3D model from 2D outlining curves in 2D slices.

In one embodiment, PerForm software may be used to provide functionalityfor loading MRI or CT scanned data, visualizing sagittal, coronal andaxial slices and scrolling over them, drawing spline curves in slices,and generating a 3D mesh model passing through a set of spline curves.In one embodiment, a 3D mesh typically is a collection of vertices,edges, and faces that may define the surface of a 3D object. The facesmay consist of triangles, quadrilaterals or other simple convexpolygons. It should be appreciated that any other curve types may beemployed instead of spline curves. For example, polyline curves may beused.

In one embodiment, a tool called “Segmentation using Landmarks” is addedto PerForm software. Such a tool provides a UI for landmarks positioningand editing. The tool also provides a button “Segment”, which invokesthe segmentation algorithm. The algorithm uses 3D image and landmarksand generates spline curves outlining the required bone.

To begin the detailed discussion of the alternative embodiment of imagesegmentation described in this section c. of the Detailed Description,wherein landmarks 777 placed on image contours are used to modify agolden bone model (e.g., golden femur or golden tibia), the resultingmodified golden bone model being segmented, reference is made to FIG.28, which is a diagram depicting types of data employed in the imagesegmentation algorithm that uses landmarks. As shown in FIG. 28, thedata employed in the segmentation algorithm 600 may be characterized asbeing two types of data. The first type of data exists in the systemonce generated and is for use with multiple patients and is notgenerated specifically for the patient for which the current imagesegmentation is being undertaken. This type of data may be called goldenmodel data 602 and is derived similar to as discussed above with respectto FIG. 11, etc. and as generally reiterated below. The golden modeldata 602 may include, for example, one or more golden femur models 603and one or more golden tibia models 604. If the joint being treated issomething other than a knee, for example, the patient's arm, then thegolden model data 602 may include another type of golden bone model, forexample, a golden radius or golden ulna.

The second type of data is specific to the patient for which the currentimage segmentation is being undertaken. This type of data may be calledinput data for segmentation algorithm 606. The input data 606 includes3D image slices data 608, which is 3D image slice data of the patientbone via MRI, CT or another type of medical imaging. The input data 606also includes landmark data 610, which is landmarks 777 positioned onboundaries of the patient bone in the image slices. The input data 606further includes patient bone characteristics 612 such as bone type(e.g., whether the bone is a tibia or femur), bone right or lefthandedness, and whether the segmentation is being done to generate anarthritic model 36 (see FIG. 1D) or a planning or restored bone model 28(see FIG. 1C). As explained below, the golden model data 602 and theinput data 606 are used in the segmentation algorithm 600 to segment the3D image employing landmarks 777.

As shown in FIG. 29, which is a flowchart illustrating the overallprocess for generating a golden femur model 603 of FIG. 28, golden femurscan image slices 616 are obtained in operation 750. For example, asdiscussed above with respect to FIG. 11 above, a representative femur618 that is free of damage and disease may be scanned via medicalimaging, such as, for example, MRI or CT. Where the golden femur model603 is to be employed in generating a bone model 22 (see block 110 ofFIG. 1C) and cartilage geometry is not of interest, the golden femurscan images slices 616 may be of a femur having damaged cartilage aslong as the bone shape is otherwise desirable (e.g., normal) and free ofdeterioration or damage. Where the golden femur model 603 is to beemployed in generating an arthritic model 36 (see block 130 of FIG. 1D)and cartilage geometry is of interest, the golden femur scan imagesslices 616 may be of a femur having both cartilage and bone shape thatare desirable (e.g., normal) and free of deterioration or damage.

The appropriate femur scan may be selected by screening multiple MRIfemur scans to locate an MRI femur scan having a femur that does nothave damaged cancellous and cortical matter (i.e., no damage in femurregions that should be present in this particular model), which has goodMRI image quality, and which has a relatively average shape, e.g., theshaft width relative to the largest part is not out of proportion (whichmay be estimated by eye-balling the images). This femur scan data,referred to herein as a golden femur scan, may be used to create agolden femur template.

It is to be appreciated that several MRI scans of a femur (or other boneof interest) may be selected, a template generated for each scan,statistics gathered on the success rate when using each template tosegment target MRI scans, and selecting the one with the highest successrate as the golden femur template.

In other embodiments, a catalog of golden models may be generated forany given feature, with distinct variants of the feature depending onvarious patient attributes, such as (but not limited to) weight, height,race, gender, age, and diagnosed disease condition. The appropriategolden mesh would then be selected for each feature based on a givenpatient's characteristics.

In operation 752 and as indicated in FIG. 30, each of the image slices616 of the representative femur 618 are segmented with a contour curveor spline 620 having control points 622 and in a manner similar to thatdiscussed above with respect to FIG. 12A, etc. For example and as shownin FIG. 30, where the golden model is to be used in the generation of abone model 22 (see block 110 of FIG. 1C) and cartilage geometry is notof interest, each segmentation region includes cancellous matter andcortical matter of the femur in a manner similar to that discussed abovewith respect to the cancellous matter 322 and cortical matter 324 of thetibia depicted in FIG. 12A, etc. Thus, as shown in FIG. 30, the contourcurve 620 excludes any cartilage matter in outlining a golden femurregion.

On the other hand, where the golden model is to be used in thegeneration of an arthritic model 36 (see block 130 of FIG. 1D) andcartilage geometry is of interest, each segmentation region the contourcurve would include cartilage matter in outlining a golden femur region.

If the golden femur scan does not contain a sufficiently long shaft ofthe femur bone (e.g., it may be desired to segment a femur in a targetMRI that may have a longer shaft), then the image segmentation can beextrapolated beyond the image to approximate a normal bone shape. Thiscan be done because the femoral shaft is quite straight and, generally,all that is needed is to continue the straight lines beyond the MRIimage, as can be understood from the extension of the contour line 620proximal of the proximal edge of the femur image 616 of FIG. 30.

In operation 754 and as illustrated in FIG. 31A, the contour curves orsplines 620 are compiled and smoothed into a golden femur mesh 624 asdiscussed above with respect to FIG. 13A, etc. As indicated in FIG. 30,in one embodiment, the segmentation curve 620 is a closed curve. Thus,the resulting golden femur mesh 624 is a closed mesh as depicted in FIG.31A.

In operation 756 and as shown in FIG. 31B, the golden femur mesh 624 isconverted into an open golden femur mesh 626, wherein the proximalportion of the golden femur mesh 624 is removed to create the opensurface model called the open golden femur mesh 626. In other words, inoperation 756 the artificial part of the femur mesh 626 is cut off,namely the proximally extending shaft portion that results from theproximal extrapolated extension of the contour line 620, so as to obtainthe open golden femur mesh 626 of FIG. 31B.

In operation 758 and as indicated in FIG. 31C, regions 628, 629 of adifferent precision are generated for the golden femur mesh 626. Forexample, when segmenting the image slices 16 for the purpose ofgenerating a golden femur mesh 626 that is used to create a 3D computergenerated bone model used to show the preoperative planning (“POP”)images to a surgeon, it is desirable that the bone geometry of the mesh626 be generated with a relatively high degree of accuracy in certainregions 628 of the mesh 626 such that the resulting 3D computergenerated bone model allows the physician to verify the POP with adesired degree of accuracy, while other regions 629 of the mesh 626 maynot be generated to such a high degree of accuracy. For example, such adegree of accuracy in the certain regions 628 of the mesh 626 can beachieved via relatively precise image segmentation. The certain regions628 of the mesh 626 having the relatively high degree of accuracy couldinclude, among others, the lower shaft area, as depicted in FIG. 31C. Inone embodiment, the relatively high accuracy of the certain regions 628of the mesh 626 should allow the physician to verify the POP within 0.5mm accuracy.

As can be understood from FIG. 31C, in one embodiment, the highprecision region(s) 628 of the mesh 626 represent a portion of thedistal anterior femoral shaft that would be contacted by the anteriorflange of a candidate femoral implant. The rest of the mesh 626 may formthe region 629 that has an accuracy that is not as precise as the highprecision region 628. Such a lower precision region 629 of the mesh 626may include the entire distal femur excluding the distal anterior regionof the shaft included within the high precision region 628. Where thegolden femur mesh 626 is employed to form other 3D computer generatedbone models, such as, for example, the bone model 22 or arthritic model36, the mesh 626 may have a different number of high precision regions628 (e.g., none, one, two, three, or more such regions 628). Also, suchregions 628 may have precisions that are greater or less than statedabove. Finally, such regions 628 may correspond to different regions ofthe bone, encompass generally the entirety of the mesh surface, orinclude other regions in addition to the region 628 depicted in FIG.31C.

While the preceding discussion regarding the open golden bone mesh isgiven in the context of the open golden bone mesh being an open goldenfemur mesh 626, as can be understood from FIG. 32A-B, the open goldenbone mesh may be an open golden tibia mesh 630 having regions 632, 633of a different precision, all of which are generated in a manner similarto that discussed with respect to FIGS. 28-31C above.

For example, as can be understood from FIGS. 32A-32B, in one embodiment,the high precision region(s) 632 of the open golden tibia mesh 630represent a portion of the proximal anterior tibial shaft immediatelydistal the tibial plateau and running medial to lateral generallyproximal the tibial tuberosity. Another high precision region 632 mayoccupy a space similar in location and size, except on the posterior ofthe tibial shaft. The rest of the mesh 630 may form the region 633 thathas an accuracy that is not as precise as the high precision region 632.Such a lower precision region 633 of the mesh 630 may include the entireproximal tibia excluding the regions of the shaft included within thehigh precision regions 632. Where the golden tibia mesh 630 is employedto form other 3D computer generated bone models, such as, for example,the bone model 22 or arthritic model 36, the mesh 630 may have adifferent number of high precision regions 632 (e.g., none, one, two,three, or more such regions 632). Also, such regions 630 may haveprecisions that are greater or less than stated above. Finally, suchregions 630 may correspond to different regions of the bone, encompassgenerally the entirety of the mesh surface, or include other regions inaddition to the regions 632 depicted in FIGS. 32A-32B.

For a discussion of an alternative embodiment of operations 250-254 ofFIG. 6, reference is first made to FIG. 33, which is a flowchartillustrating the alternative embodiment of segmenting a target bone. Inthis example, the target bone is a femur 204, but may be a tibia 210 orany other type of bone.

As indicated in FIG. 33, operation 250 obtains or, more specifically,loads the scan data (e.g., scan images 16) generated by imager 8 of thepatient's joint 14 to be analyzed. In operation 251 the landmarks arepositioned in the scan images. In other words, as can be understood fromFIG. 34, which is a flowchart illustrating the steps of operation 251,operation 251 begins with operation 251 a, wherein the images 16 arescrolled through (e.g., medial to lateral or lateral to medial) to themost medial or lateral image slice were the femur bone 204 firstappears, as shown in FIG. 35A, which, in this example, is a most lateralsagittal MRI image slice 16 where the femur bone 204 or, morespecifically, the lateral epicondyle 776 first appears. Since the slice16 of FIG. 35A is the most lateral image where bone has begun to appear,the fibula 775 can be seen adjacent the tibia 210 in such instanceswhere the image slice is positioned so as to show both the femur 204 andthe tibia 210. In operation 251 b, two or more landmarks 777 arepositioned on the outer rim of the black cortical bone 208 of the imageslice 16 depicted in FIG. 35A. As is the case with all of the imagesdepicted in FIGS. 35A-35H, in one embodiment, the landmarks are placedvia an operator sitting at a work station. In one embodiment, theoperator or user is able to add landmarks by simply clicking onto theslice image, the landmark (point) being created at the exact coordinateswhere the click has occurred. The operator is able to move existinglandmarks within the slice by selecting them and moving them with themouse, a keyboard, a pen-and-tablet system, or similar. The user is ableto delete existing landmarks by selecting them and indicating to thesoftware that they should be deleted.

In another embodiment, a touch-screen surface may be used to provideinput and display for interactive editing of landmarks and segmentationcurves. Specialized gestures may be adopted for various editingoperations.

In another embodiment, a spatial input device may be used formanipulation of landmarks, segmentation curves, and other operationsinvolving POP and jig design activities.

In operation 251 c, the image slices 16 are scrolled lateral to medialthrough approximately three slices 16 further to a new image slice 16and, at operation 251 d, it is determined if the femur bone 204 is stillvisible in the new image slice 16, which is depicted in FIG. 35B. If so,then operation 251 e adds landmarks 777 to the new image slice 16 asindicated in FIG. 35B. Specifically, as indicated in FIG. 35B, this newimage slice 16 may show the femur lateral condyle 778 and be the firstimage slice having a clear boundary 779 of the femur lateral condyle. Ascan be seen in FIG. 35B, the fibula 775 and tibia 210 are also morefully shown. Landmarks 777 are set on the clear boundary 779 of theouter rim of the dark cortical bone of the femur lateral condyle, and anadditional landmark 777 is set on the opposite side 780 on the rim ofthe black cortical bone 208. As is the case with the placement oflandmarks 777 in any of the images 16, more or fewer landmarks 777 maybe placed along the rim of the black cortical bone depicted in the image16, including landmarks being placed on the rim of the black corticalcone of the entirety of the distal femur, including the distal femurcondyle and distal femur shaft.

Operations 251 c through 251 e are repeated to set landmarks 777 at thebone contour boundaries of approximately every third image slice 16moving lateral to medial until eventually at operation 251 d it isdetermined that bone no longer appears in the present image slice. Thus,as operation 251 of FIG. 33 loops through operations 251 c-251 e of FIG.34, landmarks 777 are set at the bone contour boundaries in each of thesagittal image slices 16 depicted in FIGS. 35C-35H, which are,respectively, approximately every third sagittal image slice 16 tabbinglateral to medial through all the sagittal image slices 16 loaded inoperation 250 of FIG. 33. Thus, as shown in FIG. 35C, which represents asagittal image slice 16 approximately three slices more medial than theimage slice 16 of FIG. 35B, the femur lateral condyle 778 has a clearbone contour boundary 779, and landmarks 777 are set along the boundary779 on the rim of the dark cortical bone 208. A landmark 777 is also seton the top region 780 of the cortical bone boundary 779 where the bonecontour boundary is less clear, the landmark being positioned on the rimof the dark cortical bone 208.

As illustrated in FIG. 35D, which represents a sagittal image slice 16approximately three slices 16 more medial than the image slice 16 ofFIG. 35C, the femur shaft 781 has now appeared in an image slice 16 andboth the femur shaft 781 and femur lateral condyle 778 have clear bonecontour boundaries 779. Landmarks 777 are set along the bone contourboundaries 779 on the rim of the dark cortical bone 208.

As shown in FIG. 35E, which represents a sagittal image slice 16approximately three slices 16 more medial than the image slice 16 ofFIG. 35D, the femur lateral condyle 778 is starting to disappear, andpart of its cortical bone contour boundary 779 is not clear. Landmarks777 are only set outside the dark cortical bone 208 in the regions wherethe contour boundary 779 is clear.

As illustrated in FIG. 35F, which represents a sagittal image slice 16approximately three slices 16 more medial than the image slice 16 ofFIG. 35E, the bone contour boundary 779 has become less clear as thefemur lateral condyle 778 has decreased in size as compared to the femurlateral condyle 778 of slice 16 in FIG. 35E. The slice 16 of FIG. 30F isjust lateral of the trochlear groove 782 between the femur lateralcondyle 778 and femur medial condyle 783. The bone contour boundary 779is clear in the anterior region of the femur lateral condyle 778 and twolandmarks 777 are placed there. Additional landmarks 777 are set alongthe bone contour boundaries 779 on the rim of the dark cortical bone208.

As indicated in FIG. 35G, which represents a sagittal image slice 16approximately three slices 16 more medial than the image slice 16 ofFIG. 35F, landmarks 777 are set along the bone contour boundaries 779 onthe rim of the dark cortical bone 208. The slice 16 of FIG. 35G is inthe trochlear groove 782 between the femur lateral condyle 778 and femurmedial condyle 783. The intercondylar eminence 784 of the tibia 210 canbe seen in the slice 16 of FIG. 35G.

As indicated in FIG. 35H, which represents a sagittal image slice 16approximately three slices 16 more medial than the image slice 16 ofFIG. 35G, the femur shaft 781 has begun to disappear and the femurmedial condyle 783 has begun to appear as the slice of FIG. 35H ismedial of the trochlear groove 782 depicted in the slice of FIG. 35G.The bone contour boundary 779 is clear in the anterior region of thefemur medial condyle 783 and two landmarks 777 are placed there.Additional landmarks 777 are set along the bone contour boundaries 779on the rim of the dark cortical bone 208.

As stated above, operations 251 c through 251 e continue to be repeatedas the slices 16 continue to be tabbed through lateral to medial to setlandmarks 777 at the bone contour boundaries of approximately everythird image slice 16 until eventually at operation 251 d it isdetermined that bone no longer appears in the present image slice.Operation 251 f then scrolls medial to lateral through the image slices16 until arriving at the image slice 16 where the most medial portion ofthe femur is depicted. Operation 251 g then sets two or more landmarks777 around the bone (e.g., the medial epicondyle) in a manner similar tothat depicted in FIG. 35A with respect to the lateral epicondyle 776.This is the end of operation 251 and, as can be understood from FIG. 33,operation 252 begins by pressing the “segment” button (operation 252 a),which causes segmentation lines to be generated for each slice 16 withlandmarks 777 (operation 252 b) in a manner similar to that illustratedand discussed above with respect to FIGS. 7A-7K or as now will bediscussed below beginning with FIG. 36.

When positioning landmarks, a user needs to distribute them over thecortical bone outer surface, favoring areas where the cortical boneboundary is sharp and is more orthogonal to the slice plane,particularly favoring certain “important” areas of the bone surface(where importance is dictated by eventual contact between bone andimplant or by other requirements from POP procedure.) The user shouldonly sparsely mark up the remaining parts of the bone, particularlywhere there is a lot of volume averaging (and/or the bone surface ismore parallel to slice plane.) While the image slices depicted in FIGS.35A-35H are MRI generated image slices, in other embodiments the imagingslices may be via other medical imaging methods, such as, for example,CT.

In one embodiment, the landmark-driven segmentation algorithm describedbelow is deliberately sensitive to the number of landmarks (points)placed at a given area of the bone. So for instance, if the user desiresthe auto-generated bone mesh to very accurately pass through particularspots on the slice, the user can place more than one landmark on thatsame spot or very near that spot. When there is a high concentration oflandmarks in a small area of the bone, the auto-generated mesh will bebiased to more accurately model that area. The software indicates to theuser, making it visible at a glance whenever more than one landmark islocated within the same small area of the image.

In one embodiment, instead of putting landmarks in every three slices, auser may position landmarks in every slice but use three times fewerlandmarks in each slice. The result of the segmentation usually variesvery little depending on how a user distributes landmarks around thebone surface as long as the entire surface is covered.

While much of the following discussion takes place in reference to thesegmentation of a femur (operation 252 of FIG. 6), the conceptsdiscussed herein are readily applicable to the segmentation of a tibia(operation 258 of FIG. 6). Additionally, the concepts discussed hereinare readily applicable to both the left or right knee. Different goldentemplate data may be used to segment the left tibia, right tibia, leftfemur or right femur for bone models 22 or planning models 28.Additionally, other embodiments may segment other models and or joints,including but not limited to, arthritic models 36, hip joints, elbowjoints, etc. by using an appropriate golden template of the feature ofinterest to be segmented.

As shown in FIG. 36, which is a flowchart illustrating the process ofsegmenting the target images 16 that were provided with landmarks 777 inoperation 251, the full or entire golden femur mesh 626, including itsregions 628, 629 in FIG. 31C, is deformed in operation 770 to matchlandmarks 777 and appropriate features, such as, for example, the outeredges of dark cortical bone, in the target scan images 16.

As discussed below with respect to FIG. 37, a method is provided formapping the golden femur mesh into the target scan using registrationtechniques. Registration may be thought of as an optimization problemwith a goal of finding a spatial mapping that aligns a fixed object witha target object. Generally, several registration operations may beperformed, first starting with a low-dimensional transformation group tofind a rough approximation of the actual femur location and shape in thetarget image. This may be done to reduce the chance of finding wrongfeatures instead of the femur of interest. For example, if a free-formdeformation registration was initially used to register the golden femurmesh to the target scan data, the template might be registered to thewrong feature, e.g., to a tibia rather than the femur of interest. Acoarse registration may also be performed in less time than a fineregistration, thereby reducing the overall time required to perform theregistration. Once the femur has been approximately located using acoarse registration, finer registration operations may be performed tomore accurately determine the femur location and shape. By using thefemur approximation determined by the prior registration operation asthe initial approximation of the femur in the next registrationoperation, the next registration operation may find a solution in lesstime. It is to be understood that similar considerations apply tosegmentation of other entities (and not just the femur.)

In one embodiment, each registration operation may employ a registrationframework. The registration framework may be based on three generalblocks. The first block defines a transformation model (or a class oftransforms) T(X), which may be applied to coordinates of a fixed (orreference) object (e.g., a golden femur template) to locate theircorresponding coordinates in a target image space (e.g., an MRI scan).The second block defines a metric, which quantifies the degree ofcorrespondence or similarity between features of a fixed (or reference)object and the target object (that is landmarks and appropriate targetimage features) achieved by a given transformation. It should be notedthat instead of a metric that defines the degree of correspondence, anopposite to it function is defined, which is call the defect function.The third block defines an optimization algorithm (optimizer), whichtries to maximize the reference and the target objects similarity (orminimize the opposite defect function) by changing the parameters of thetransformation model. Thus, as discussed below in detail with referenceto FIG. 37, in every registration operation 770 a-770 c and 770 e thereis a need to specify three blocks: (1) class of transforms; (2) metric(or defect) function; and (3) optimization algorithm. In one embodiment,the same third block may be used in all four registration steps. Forinstance, a gradient descent optimizer or conjugate gradient descendoptimizer may be used. Alternatively, any other appropriate optimizationalgorithm, such as Monte Carlo, simulated annealing, genetic algorithms,neural networks, and so on, may be used.

As shown in FIG. 37, which is a flowchart illustrating the process ofoperation 770 of FIG. 36, in operation 770 a translation transforms areused to register the full or entire open golden femur mesh 626 to thelandmarks 777. More specifically, in operation 770 a, the open goldenfemur mesh 626 may be approximately registered to landmarks 777 using acoarse registration transformation. In one embodiment, this may be doneby finding appropriate translation transform parameters that minimizetranslation misalignment with landmarks of the reference open goldenfemur mesh mapped onto the target femur of the target image, wherelandmarks 777 are positioned. This coarse registration operationtypically determines an approximate femur position in the MRI scan.During this operation, the reference open golden femur mesh 626 may beoverlapped with the target femur of the target image using a translationtransformation to minimize translational misalignment of the femurs. Atranslation transform, translates (or shifts) all the points in 3D spaceby the same 3D vector. That is, the reference femur may be mapped intothe target image space by shifting the reference open golden femur meshalong one or more axes in the target image space to minimizemisalignment. During this operation the reference object is not rotated,scaled or deformed. In one embodiment, three parameters for thetranslation transformation may be generated: one parameter for eachdimension that specifies the translation for that dimension. In oneembodiment, the final parameters of the translation transform minimizingthe misalignment of the mapped reference femur image coordinates intothe target image space may be found using a gradient descent optimizer.In other embodiments, other types of optimizers may be utilized, such asfor instance an Iterative Closest Point (ICP) algorithm.

Optimization of mesh alignment with respect to landmarks is based onminimizing a cost function D, which in one embodiment can be the sum,across all landmarks, of the squared distance from each landmark point777 to the transformed open golden mesh. The same cost function may beused for steps 770 a-770 c. Methods for computing this cost function andits gradient are covered in more detail later in this DetailedDescription.

After an optimal transform has been found, it is applied to all thegolden femur data (i.e., the closed golden femur mesh 624, open goldenfemur mesh 626, and golden femur mesh regions 628, 629. The nextoperation (i.e., operation 770 b of FIG. 37, which is discussedimmediately below) is then started with transformed golden femur data.As can be understood from the following discussion, after everyconsecutive operation 770 a, 770 b, 770 c and 770 e of FIG. 37, thetransform found during the registration step is applied to all thegolden femur data. As a result, after each operation the golden femurdata is successively made more closely aligned with the femur in thetarget image.

In operation 770 b of FIG. 37 similarity transforms are used to registerthe full or entire open golden femur mesh 626 to the landmarks 777.Specifically, operation 770 b further refines the object's registrationdetermined by operation 770 a. This may be done by approximatelyregistering the open golden femur mesh 626 to landmarks 777 using asimilarity transformation. In one embodiment, a similaritytransformation may be performed in 3D space. The reference open goldenfemur mesh may be rotated in 3D, translated in 3D and homogeneouslyscaled to map its coordinates into the target MRI scan data to minimizemisalignment between the open golden femur mesh and the landmarks in thetarget MRI scan. In some embodiments, a center of rotation may bespecified so that both the rotation and scaling operations are performedwith respect to the specified center of rotation. In one embodiment, a3D similarity transformation, specified by seven parameters, may beused. One parameter specifies the scaling factor, three parametersspecify a versor that represents the 3D rotation, and three parametersspecify a vector that represents the 3D translation in each dimension. Aversor is a unit quaternion that provides a convenient mathematicalnotation for representing rotations of objects in three dimensions.

In one embodiment, local minimization techniques may be employed withthe similarity transformation to obtain a refined registration of thereference open golden femur mesh onto the target MRI scan that is notfar from the registration of the reference open golden femur mesh ontothe target MRI scan found in previous operation 770 a and used as theinitial starting approximation. For instance, gradient descent,conjugate gradient descent, or ICP optimization may be used. After thebest transform is found for operation 770 b of FIG. 37, the transform isapplied to the golden femur data in a manner similar to that ofoperation 770 a.

In operation 770 c of FIG. 37 affine transforms are used to register thefull or entire open golden femur mesh 626 to the landmarks 777.Specifically, operation 770 c further refines the image registrationdetermined by operation 770 b. In one embodiment, an affinetransformation may be used to register the open golden femur mesh 626 tolandmarks 777 in the target MRI scan data. In one embodiment, theapproximate femur registration found during operation 770 b may be usedas the initial starting approximation for the affine transformation ofoperation 770 c.

An affine transformation typically is a linear transformation followedby a translation. The affine transformation preserves collinearitybetween points (i.e., three points which lie on a line continue to becollinear after the transformation) and ratios of distances along aline. In one embodiment, a 3D affine transformation, specified by 12parameters, may be utilized. Nine parameters of the affinetransformation specify the linear transformation (which may berepresented by a three by three matrix) and three parameters of theaffine transformation specify the 3D translation in each dimension. Theparameters of the affine transform that minimizes the misalignment ofthe open golden femur mesh with landmarks may be found using again localminimization techniques, such as gradient descent or conjugate gradientdescent optimization.

After the best transform is found for operation 770 c of FIG. 37, thetransform is applied to the golden femur data. The transformed goldenfemur data from operation 770 c is then employed in the preparatory stepof detecting appropriate image edges, namely, operation 770 d, which isdiscussed below. Those edges will be later used in operation 770 e ofFIG. 37, as discussed below. The transformed golden femur data fromoperation 770 c is used as reference data similar to the previousoperations.

A discussion of image edges is now provided before discussing thedetails of operation 770 d of FIG. 37. Image edges consist of thosepoints in the image where the image contrast significantly changesbetween neighbor pixels (or voxels) and this contrast change isconsistent along several neighboring points distributed over a smoothcurve. For example, points that lie between the light cancellous bonepixels and dark cortical bone pixels form an image edge. Similarly, thepoints that lie between the dark cortical bone pixels and the grayishcartilage pixels form an image edge. Yet a configuration involving aone-pixel black spot and the surrounding light pixels does not form animage edge because the light points represent a curve with too muchcurvature, whereas the dark point represents a curve that is toodiscontinuous (spanning only a single voxel.) Usually there is an imageedge that separates one type of the body tissue from a neighboringdifferent type of body tissue.

The purpose of segmenting an image is to be able to separate in theimage certain body tissues from the surrounding tissues. Ideally, thesegmentation boundaries (or curves) should lie mostly in the imageedges. A general MRI or CT image contains lots of edges separatingvarious body tissues from the neighboring tissues. Yet when segmenting,there is only interest in certain tissues and thus particular edgesonly. Operation 770 d is intended to find those edges that are ofinterest for segmenting a particular body object.

In particular in case of the segmentation of any of the versions of thefemur planning model 22, 28 (shown in blocks 110 and 115, respectively,of FIG. 1C), operation 770 d of FIG. 37 will find the edges thatseparate the cortical femur bone from the outside knee tissues (i.e.,the tendons, ligaments, cartilage, fluid, etc.). In some embodiments,operation 770 d will not find the edges that separate the femurcancellous bone from the femur cortical bone. In other embodiments,operation 770 d will find the edges that separate the femur cancellousbone form the cortical bone.

Operation 770 d may also find some edges that are of no interest to thefemur planning segmentation. Most of those edges of no interest will lieat significant distance from the femur boundary surface and, as aresult, the edges of no interest will not influence the next operationin the algorithm, namely, operation 770 e of FIG. 37.

In some cases, some of the edges of no interest might happen to be veryclose to the edges of interest. Such nearby edges of no interest arelikely to be the edges separating the cartilage tissue from the othertissues outside the bone. Such edges might adversely influence the nextoperation in the algorithm, namely, operation 770 e of FIG. 37, and leadto inaccurate segmentation. In some embodiments, this inaccuracy can beremedied by the user providing extra landmarks 777 in the area that islikely to cause such inaccuracies or manually fixing the spline curvesduring the verification and adjustment operations.

The result of the operation 770 e of FIG. 37 will be a 3D image of thesame size as the target scan data. The resulting 3D image can be calledan edges image. The voxels in the edges image correspondent to strongedges will have highest intensities, the non-edge voxels will have lowintensities, and the voxels correspondent to weak edges will haveintermediate intensities. Discussion of the operation 770 d of FIG. 37is now provided.

In operation 770 d of FIG. 37 appropriate edges of the target images aredetected near the transformed open golden femur mesh 626. For example,as indicated in FIG. 38A, which is a flowchart illustrating the processof operation 770 d of FIG. 37, in operation 770 d 1 the signed distanceimage is computed for the transformed golden femur mesh 626. A signeddistance map is a distance map of a region in 2D (or 3D) and is afunction in 2D (or 3D). The signed distance value for a point equals thedistance from the point to the boundary of a region. A signed distancevalue can have a positive or negative value. For example, when a pointis inside the region, the signed distance value of the point is thedistance from the point to the boundary of the region in the form of anegative value. When a point is outside the region, the signed distancevalue of the point is the distance from the point to the boundary of theregion in the form of a positive value. If the signed distance mapfunction is computed in a regular grid of points in 2D (or 3D)correspondent to image pixels (or voxels) and stored as a 2D (or 3D)image representation, the result can be said to be a 2D (or 3D) signeddistance image.

Thus, from the preceding discussion, it can be understood that thesigned distance for a watertight surface is a function that has absolutevalues equal to the regular (Euclidean) distance, but the values alsohave a sign. The sign is negative for the points inside the surface, andthe sign is positive for the points outside the surface. The open goldenfemur mesh 626 transformed in operations 770 a-770 c of FIG. 37 is usedin operation 770 d 1 of FIG. 38A. By the time of operation 770 d 1, theopen golden femur mesh 626 may quite closely match the landmarks 777positioned in the target image and, as a result, the open golden femurmesh 626 also matches quite closely the target femur bone in the targetimage. Since the golden femur mesh 626 is a watertight mesh, the maskimage marking may be computed as “1” for all voxels that lie inside theopen golden femur mesh 626 and as “0” for all the voxels that lieoutside the mesh. The Signed Danielsson Distance Map Image Filter fromthe ITK library can then be used to compute the signed distance to themask boundary, which is approximately the same as the signed distance tothe mesh. It may be desired to have greater accuracy close to the mesh.If so, then for the voxels where the absolute value of the signeddistance is small, the distance to the mesh may be recomputed by findingthe closest points via a more exact method, as detailed later in thisspecification.

In operation 770 d 2 the gradient of the signed distance image iscomputed. As can be understood from FIG. 38B, the gradient of the signeddistance image contains a vector 1000 in every voxel. The vector 1000represents the gradient of the signed distance image at the particularpoint of the voxel. Because the signed distance image represents thesigned distance to the transformed open golden femur mesh 626, whichfollows closely the boundary of the femur bone in the target image, thegradient image has gradient vectors nearly orthogonal to the boundary ofthe target femur in the voxels close to the boundary.

The contour line 626 in FIG. 38B represents the approximate segmentationmesh surface found in the previous registration step of operation 770 cof FIG. 37. The vectors 1000 show the gradient of the signed distancefor the contour 626. The starting end of the vector 1000 is the point orvoxel where the vector 1000 is computed. The gradient of a signeddistance has a vector direction in every point or voxel toward theclosest point in the contour 626. Vectors 1000 are oriented from insideto outside the contour 626. Each vector 1000 has a unit length.

In operation 770 d 3 the gradient of the target image is computed. Ascan be understood from FIG. 38C, which is an enlarged view of the areain FIG. 38B enclosed by the square 1002, the gradient of the targetimage has gradient vectors 1004 orthogonal to the edges 1006, 1008 inthe target image, and the length of those vectors 1004 is larger forstronger edges and smaller for weaker edges. Such vectors 1004 arealways oriented from the darker image region to the lighter image regionor, in other words, from darker pixels towards brighter pixels. Thevectors 1004 are longer where the contrast is higher. For purposes ofillustration in FIG. 38C, the vectors 1004 illustrated are only longvectors corresponding to high contrast pixels associated with strongedges. The gradient vectors 1004 can be used to identify the outercortical bone boundary 1006 and other edges 1008, 1010 that are not ofinterest for the analysis.

Finally, operation 770 d of FIG. 37 is completed via operation 770 d 4of FIG. 38A, wherein the edges image is computed by correcting thegradient of the target image with the gradient of the signed distanceimage. As can be understood from FIG. 38D, the edges image is computedby combining the results from operations 770 d 2 and 770 d 3. Dependingon the type of 3D computer generated bone model being generated from thesegmented images, different boundary edges may be of relevance. Forexample, if the images are being segmented to generate a bone model 22,the boundary edges that are of interest contain dark cortical voxelsinside and lighter cartilage or other voxels outside. As a result, thevoxels that are of interest are those voxels that have similarlyoriented gradients 1000, 1004 computed in operations 770 d 2 and 770 d 3as shown in FIGS. 38B and 38C, respectively. In every voxel the vector1004 from operation 770 d 3 is projected onto the vector 1000 fromoperation 770 d 2. When the projection of image gradient vector onto asigned distance gradient vector points in the same direction as thesigned distance vector, its magnitude is taken as the voxel value forthe resulting edges image. When it points in the opposite direction (orhas no magnitude at all), “0” is taken as the voxel value for theresulting edges image.

The resulting edges image has high values in the target femur corticalbone outer boundary 1006. However, the edges image does not have manyother high values close to the transformed open golden femur mesh withone exception, namely, the voxels on the boundary between the targetfemur cartilage and the outsight bright voxels (for example fluidvoxels) might have high values.

As can be understood from FIG. 38D, the gradient of the signed distancevectors 1000 are uniformly oriented orthogonal to the bone surface andgo from inside to outside of the bone. The image gradient vectors 1004are oriented differently near different image details. For the points inthe bone outer boundary 1006, the vectors 1004 are almost parallel tothe vectors 1000. For two of the other boundaries 1008, the vectors1000, 1004 are generally oppositely oriented. Finally, for the lastboundary 1010, the vectors 1000, 1004 are quite differently orientedfrom each other. As a result of the combination of the vectors 1000,1004 and an analysis of their directional relationships to each other,the edges image will be for FIG. 38D as follows. For points associatedwith the bone contour line 1006, the edges image will reflect the lengthof the image gradient vector. For points associated with contour lines1008, the edges image will be zero. For points associated with thecontour line 1010, the edges image values will be smaller than thelength of the image gradient vector associated with the bone contourline 1006. Thus, the edges image will tend to have the largest valuesfor the points of the bone contour line 1006.

The high values correspondent to a cartilage/fluid boundary mightnegatively impact operation 770 e of the registration in FIG. 37.Consequently, it may be desirable to suppress those values. This can bedone in the beginning of operation 770 d 3 of FIG. 38A. Specifically, awindowing filter may be applied to the whole target image. A window [w0,w1] may be used, where w0 will be the minimum value in the image, and w1will be approximately the value correspondent to the cartilageintensity. The filter will replace the high intensity values in theimage with w1 value, and thus the boundary between the cartilage and thelighter matters will disappear. For the type of MRI images that may beused, the w1 value correspondent to the median of all the values in theimage works quite well. Although such a filter may not always suppressthe cartilage boundary entirely, it makes cartilage outer boundary verymuch weaker in the image and, as a result, the cartilage has less of animpact in operation 770 e of FIG. 37.

In one embodiment, a more sophisticated method for suppressing thecartilage boundary may be employed. The cartilage intensity values maybe estimated by comparing the voxel values near landmarks 777 along thesigned distance gradient direction. The values before a landmarkcorrespond to the cortical bone intensities, while the values after thelandmark correspond to the cartilage intensity. Thus for every landmark,a value may be found that represents an “Out of cortical bone”intensity. Such values may be interpolated into the whole image and thiswindowing function may be applied rather than the constant windowingvalue w1.

It should be appreciated that a lesser resolution than the target imageresolution may be used in all the images participating in the edgesimage computation. For example, an in-slice voxel size of 1 mm may beused rather than −0.3 mm in the target image. Using coarser resolutionin effect smoothes out the data, allowing a more stable edgescomputation. It also significantly speeds up the computation. In case ofvery noisy target images, an additional smoothing step may be applied.

Operation 770 of FIG. 36 is completed via operation 770 e of FIG. 37,wherein the full or entire golden femur mesh 626, including its regions628, 629, are simultaneously registered to landmarks 777 and image edgesrespectively using B-spline deformable transforms. Specifically,operation 770 e of FIG. 37 further refines the image registration of theboundary golden femur region. In one embodiment, a spline transformationmay be used to register the open golden femur mesh 626 into the MRI scandata (target image space). In one embodiment, 3D B-Spline deformabletransforms may be employed.

A B-Spline deformable transformation typically is a free formdeformation of an object using a deformation field where a deformationvector is assigned to every point in space. For example, a 3D B-splinedeformable transform T may specify a 3D vector V(P) for every point P inthe original 3D space that is moved by T such that T:P→P+V(P).

In one embodiment, a B-Spline transformation may be specified with M×Nparameters, where M is the number of nodes in the B-Spline grid and N isthe dimension of the space. In one embodiment, a 3D B-Spline deformabletransformation of order three may be used to map every reference image3D point into the target MRI scan by a different 3D vector. The field ofvectors may be modeled using B-splines. Typically a grid J×K×L ofcontrol points may be specified where J, K, and L are parameters of thetransformation.

In one embodiment, splines of order three may be used with a grid27×9×11 of control points. That is, the transformation employs 27control points in the medial/lateral direction (i.e., the x direction),9 control points in posterior/anterior direction, and 11 control pointsin distal/proximal direction. Two control points in each dimension(i.e., 2 of 27 in the x direction, 2 of 9 in the y direction and 2 of 11in the z direction) may be used to specify boundary conditions. As such,the inner spline nodes may form a grid of size 25 by 7 by 9 and theboundary conditions increase the grid to size 27 by 9 by 11. Theparametric set for this transformation has a dimension of 3×27×9×11=8019(i.e., at each node of a 27×9×11 grid of control points, there isspecified a 3-dimensional transformation vector; a nonlinearinterpolation of transformation vectors for points located betweenadjacent nodes, is governed by spline equations.) The final parametersof the spline transformation that minimizes the misalignment between thereference golden femur template and the target MRI scan data may befound.

In operation 770 e of FIG. 37 a different metric (or defect function)may be used as compared to what was used in operations 770 a, 770 b, and770 c. Specifically, a combined defect function may be used. Thecombined defect function may be defined as a linear combination of thedefect function D (same as in operations 770 a, 770 b, and 770 c) anddefect functions D_i that evaluate the discrepancy between the goldenmesh regions 628, 629 and the scan image edges defined in operation 770d of FIG. 37.

The defect function D_i, or rather its opposite metric functionM_i=−D_i, for a given Golden Mesh Region R_i may be defined as follows.All the vertices in the golden mesh region R_i, are taken, a transformis applied to them, and the correspondent intensities are evaluated inthe edges image. M_i may be set to be the sum of those intensities.Thus, when more vertices from the transformed golden mesh region R_icome close to the image edges, a higher metric value is the result.

When defining the combined metric or defect, that is when taking thelinear combination of D and all the D_i, the coefficients in the linearcombination need to be specified. It may be desirable to use a very highcoefficient with D because we want to follow very precisely thelandmarks 777 provided by a user. Smaller coefficients may be employedwith D_i. The latter coefficients might be also different. The highercoefficients may be used for those regions of the bone that require agreater degree of precision, the associated image segmentation needingto result in more clearly defined regions. The lower coefficients may beused for those regions of the bone that do not require a high degree ofprecision, the associated image segmentation resulting in less clearlydefined regions.

Some bones may have a higher degree of shape variations across thepopulation than is found with the knee region of the femur. For example,the shape of the tibia may vary more from patient to patient than doesthe shape of the femur. As a result, the affine transformation may notprovide a close enough registration of the golden tibia template to thetarget tibia in the target scan. This may cause the Splinetransformation to find a local optimum that may be far from the actualtibia in some areas. In one embodiment, an additional registrationoperation between the affine transform and spline transform operationsmay be performed to more closely align the golden tibia and the targettibia, allowing the spline transform to converge to the correct localoptimum rather than a nearby, but wrong, local optimum.

The class of transforms utilized generally should allow more flexibility(or degrees of freedom) than the Affine transform and less flexibilitythan the B-spline transforms. The number of degrees of freedom generallyis equal to the number of transform parameters. In one embodiment, aclass of transforms with more than 12 parameters and less than 3×27×9×11parameters may be used. For example, a B-spline transform with fewercontrol points than used in the subsequent spline transform may be usedfor the additional transform operation. Alternatively, the deformationsmay be modeled using quadric rather than cubic functions.

In another embodiment, several golden tibia templates may be used thatrepresent typical tibia variations, e.g., golden tibia templates forvarus, valgus, and normal tibia. In one embodiment, each of the goldentibia templates may be used during the translation, similarity andaffine transform registration operations to find the template thatprovides the best match (e.g., best correlation) in the affine transformregistration operation. This template may then be used in the remainingregistration operations.

Finally, in one embodiment, the tibia registration may be improved byperforming the tibia segmentation after the femur segmentation andadding a restriction on the tibia registration transformations such thatthe tibia may not penetrate the femur. In one embodiment, this may beimplemented by introducing a penalty for the penetration. In the targetMRI all the voxels that lie inside the femur segmentation curves may bemarked. The metric functions, described in more detail below, that areused in the registration operations may be modified to include a penaltyterm. The penalty term may be computed by taking points in the goldentibia mesh, applying a transform to the set of points, determining if atransformed sample point falls into any of the marked voxels, and addinga large value to the penalty term for each transformed sample point thatfalls into any of the marked voxels.

In each of the above registration operations, a metric may be used toquantify the degree of correspondence between the reference objects andtarget image achieved by a given transformation. In one embodiment, themetric quantitatively measures how well the transformed golden femurdata fits the target image (e.g., a target MRI scan) and landmarkspositioned there.

As discussed above, metrics M=−D, M_i=−D_i, and their linear combinationare used in operations 770 a-770 d of the registration. An explanationis now given regarding the details on how to compute those metrics. Asfar as using those metrics with optimizers that require computations ofthe gradient of the metric, it is also explained how to compute thegradient of those metrics.

When computing the metric M or rather the defect D, such a computationcan include finding the sum of the squared distances from each landmarkpoint 777 to the transformed open golden mesh. In order to make thiscomputation as quickly and efficiently as possible, the following can bedone. First, a B-Spline transformation of a mesh is no longer a mesh.The plane triangles forming the original mesh get curved over thetransformation, and the triangles are no longer planar. Rather thancomputing distances to curved triangles, which would be verycomputationally expensive, planar triangles connecting the transformedvertices are used. Very little precision is lost with this substitutionbecause the triangles are very small.

Next, after finding the transformed mesh, it is desirable for everyLandmark point to find the closest point in the transformed meshtriangles and take the squared distance to it. A spatial subdivisionscheme is used to sort all the triangles by spatial location. An octreesubdivision is used, although other schemes (kd-tree, fixed size grid,etc.) would work as well. The spatial subdivision helps to find aclosest mesh triangle and a closest point in it using an order of LOG(n)operations where n is the number of triangles in the mesh.

The optimizers used in the registration steps require the computation ofthe gradient of the metric function, which depends on the appliedtransform, over the transform parameters.

In one embodiment the metric function may be a composition of severalfunctions. For the metric M (or cost function D), for example, thefollowing functions are used in the composition: a) mesh transformation,b) distance from a Landmark point to the transformed mesh, c) squareddistance, d) sum of squares, e) inverse of the sum.

For a composition of functions, determining the gradient involvesfinding partial derivatives for each function and then applying thechain rule. The derivatives of the algebraic functions are computed bystandard formulae. The only non-trivial computation in the abovefunctions is the computation of the partial derivative of a distancefrom a point to the transformed mesh.

For the latter computation, it may involve using an approximate method.Namely, take the closest triangle found in the metric computationfollowed by taking the plane containing that triangle. This planeapproximates the transformed mesh surface in some small neighborhood ofthe closest point. One of the transform parameters is changed by a verysmall amount. It is observed where the former closest triangle is mappedafter the variation.

The plane containing the varied triangle is taken. This planeapproximates the varied transformed mesh surface in some smallneighborhood of the varied closest point. The distance from the landmarkpoint to this varied plane is taken. It is approximately the distancefrom the landmark point to the whole varied transformed mesh surface.Now the difference between the varied distance and the original distanceis taken and divided by the value of the parameters variation. Thisgives approximately the partial derivative for this parameter.

In order to compute the gradient of the metric D_i with respect to thetransform parameters, the gradient image of the edges image is computedright after the computation of the edges image itself. To compute thepartial derivative of D_i over a transform parameter, the computationmay take place for the derivative of every transformed vertex motionover that parameter using the chain rule. This derivative will be a 3Dvector. Its dot product is taken with the correspondent Gradient vectorof the Gradient Image of the Edges Image and the values are summed allover the vertices.

Finally, since the combined defect function is a linear combination ofdefect functions D and D_i, then the gradient of the combined defectfunction with respect to a given transform, is correspondingly a linearcombination (with the same coefficients) of the gradients of D and D_iwith respect to that same transform.

In summary and as can be understood from the immediately precedingdiscussion regarding operations 770 a-770 c and 770 e of FIG. 37,translation transforms (operation 770 a), similarity transforms(operation 770 b), affine transforms (operation 770 c), and B-Splinedeformable transforms (operation 770 e) are employed as part ofaccomplishing operation 770 of FIG. 36. Because in operations 770 a-770c of FIG. 37 it is intended to register the open golden femur mesh tolandmarks, the metric (or defect) function should evaluate how close thetransformed open golden femur mesh is to landmarks. Thus, in operations770 a-770 c, there is a selection of the defect function D to be the sumof squared distances from landmarks to the deformed open golden mesh. Inthe operation 770 e, a simultaneous registering of several parametersmay be defined as a combined metric that will take into account all theparameters. The combined defect function may be defined as a linearcombination of the defect function D (same as in operations 770 a-770 c)and defect functions D_i that evaluate the discrepancy between thegolden mesh regions and the scan image edges defined in operation 770 dof FIG. 37.

Once operation 770 of FIG. 36 is completed, the process of FIG. 36 thencontinues with operation 772, wherein the deformed open golden femurmesh 626 and associated regions 628, 629 are segmented followed byoperation 773, wherein the resulting segmentation curves areapproximated with splines. The process continues with operation 774,wherein the contour lines or splines generated in operation 773 aremodified to match landmarks.

In other words, in operation 774 the segmentation curves are refined tomore precisely match the landmarks. Typically, the segmentation curvescreated via the above-described algorithms match the landmarks quiteclosely. Accordingly, most any simple algorithm for a local curveadjustment can work to further refine the precision of the match of thesegmentation curves with the landmarks.

In one embodiment of operation 774 when further refining thesegmentation curves to match landmarks, only those curves that belong toslices that contain landmarks are considered. When a curve belongs to aslice with landmarks, it is assumed that it should rather precisely gothrough all the Landmarks. In one embodiment, a curve may be consideredto be precisely enough located relative to a landmark if its distancefrom the landmark (“Tol”) is Tol=0.3 mm or less. Most often all thelandmarks are within the Tol distance from the curve. However, sometimesa few of the landmarks are further than the Tol distance from the curve.As can be understood from the following discussion regarding operation774, for every curve generated via the above-described algorithms, eachlandmark in a slice is iterated. If a landmark is not within Toldistance from the curve, a correction algorithm is applied to the curveas described below with respect to operation 774.

Operation 774 locally modifies the spline curve to fit a selectedlandmark. Specifically, as can be understood from FIG. 39, which is aflowchart illustrating the process of operation 774 of FIG. 36 in thecontext of a particular image slice containing landmarks and asegmentation contour, in operation 774 a a landmark 777 is identified.In operation 774 b, the distance of the identified landmark 777 to thespline generated from the golden femur mesh 626 is computed. Morespecifically, the algorithm of operation 774 first computes distancesfor all the other landmarks in the same slice to avoid making thedistance relationships of the landmarks and curve worse.

In operation 774 c, an arc of the contour line or spline that is theclosest to the landmark is identified. Specifically, the closest splinearc [A, B] to the selected landmark L is located, where A and B areconsecutive vertices in the spline curve.

In operation 774 d, the arc is modified to include the identifiedlandmark, resulting in a modified contour line or spline. Specifically,the vertices A and B are moved iteratively so that the arc [A, B] fitsL. For each iteration, the closest point C in [A, B] to L is found. Theratio α: (1−α) is then found in which C divides [A, B]. Next, A is movedby (1−α)*(L−C), and B is moved by α*(L−C). The process stops when0.5*Tol distance is achieved.

In operation 774 e, distances of other landmarks to the modified splineare computed and reviewed in operation 774 f to verify operations 774a-774 d have not made the fit between the spline and other landmarksworse. In other words, the following is checked for every landmark.First, it is checked to see if the new distance between the spline andlandmark is within Tol and, if it is, then the relationship between thespline and landmark is acceptable. Second, it is checked to see if thenew distance between the spline and landmark is smaller than the olddistance and, if it is, then the relationship between the spline andlandmark is acceptable. Third, it is checked to see if the new distanceis higher than Tol, the old distance was higher than Tol, and the newdistance increased by less than 0.5*Tol. If the answer is yes withrespect to all three of the elements of the third check, then therelationship between the spline and landmark is acceptable. For all theother cases, the relationship between the spline and landmark is notacceptable. If the distance relationships between the spline and all ofthe landmarks are considered acceptable, the process outlined in FIG. 39is completed for the identified landmark, and the process of FIG. 39 canthen be run for another identified landmark until all landmarks havegone through the process of FIG. 39.

If any of the distance relationships between any landmark and the splineare found to be unacceptable in operation 774 f due to a modification ofthe spline with respect to a selected landmark according to operations774 a-774 d, then per operation 774 g the spline modification fromoperation 774 d is disregarded and a more local modification isemployed, as described below with respect to operations 774 h-774 k ofFIG. 39. The more local modification will add a new vertex into thespline making the spline more flexible in this region. The more localmodification will then move the new vertex to L, and this will affect avery small area of the spline. Thus, the chance of decreasing the fit toother landmarks will be very small. The more local modification occursas follows.

In operation 774 h, a point in the identified arc closest to theidentified landmark is identified. Specifically, the point C in arc [A,B] that is the closest to landmark L is found.

In operation 774 i, a spline vertex is inserted at the identified pointC. With the insertion of a new vertex, the spline shape usually changesin the two immediately adjacent neighbor arcs on both sides of the arc[A, B]. As a result, the arc spline can become too wavy in the vicinityof the arc [A, B].

To remedy the situation, the arc [A, B] is adjusted to fit the originalspline in operation 774 j. Specifically, the vertices A and B aremodified to try to fit the new spline as closely as possible to theoriginal spline. In doing so, a measure of closeness (i.e., how closelythe new spline follows the original spline in the six neighboringarcs—three to each side of the new control point C) may be computed asfollows. In one embodiment, the six spline arcs are sampled such thatthere are twenty or so sample points in every arc of the spline (i.e.,20*6 sample points). Then, the sum of the squared distances from thesample points to the original spline may be computed. Next, thecoordinates of th\e A and B vertices (control points) are varied (i.e.,two parameters for each of A and B, that is four parameters). Then, alocal optimization algorithm is used to find the closest spline. Thisprocess may be similar to the process of fitting a spline to a polyline,as described elsewhere in this Detailed Description.

Per operation 774 k, the identified point is moved to the identifiedlandmark. Specifically, the spline vertex C is moved into the landmarkpoint L.

The process outlined in FIG. 39 is completed for the identifiedlandmark, and the process of FIG. 39 can then be run for anotheridentified landmark until all landmarks have gone through the process ofFIG. 39.

Once the process of FIG. 39 is completed for all landmarks and theassociated contour lines or splines, the process of operation 774 ofFIG. 36 is considered complete, which completes the process of FIG. 36for the operation 252 a of FIG. 33. The process of operation 252 a inFIG. 33 is now complete. The image slices 16 are then scrolled over toverify if the segmentation results are acceptable, as indicated byoperation 252 c. In operation 253, if the segmentation is acceptable,then the segmentation process of FIG. 33 ends.

As can be understood from FIG. 33, if in operation 253 the segmentationis not acceptable, then the segmentation of each offending slice 16 ismodified by adding additional landmarks 777 and/or modifying thelocations of existing landmarks 777 per operation 254 of FIG. 33. Forexample and as can be understood from FIG. 40, a first spline 800, whichis generated via a first run through operation 252 of FIG. 33, hascontrol points 802 and extends along first landmarks 777 a placed in theslice 16 of FIG. 40 during operation 251 of FIG. 33.

During operation 253 of FIG. 33 the segmentation of the slice 16 of FIG.40 is identified as being unsatisfactory in the location called out byarrow A in FIG. 40. A new landmark 777 b is added in the location calledout by arrow A per operation 254 and operation 252, or morespecifically, operations 774 b-774 e of the algorithm of FIG. 39, arerepeated to generate a second spline 804, which has control points 806and extends along both the first landmarks 777 a and the second landmark777 b. As can be understood from FIG. 40, the first spline 800 and thesecond spline 804 are generally identical and coextensive, except in theregion identified by arrow A. The segmentation of the second spline 804is then approved or disapproved per operation 253. If approved, then thesegmentation process of FIG. 33 ends. If disapproved, then the secondspline 804 is further modified per operation 254 in a manner similar toas discussed above with respect to FIG. 40.

In one embodiment of operation 254 of FIG. 33, the spline may besimultaneously modified near a new added landmark or near movinglandmarks to fit the moving landmarks. In doing so, it may be the casethat the user is satisfied with the corrected splines. As a result, theprocess of FIG. 33 may simply end at operation 254 as if the entirety ofoperation 252 had been completed and the segmentation was foundacceptable at operation 253.

In one embodiment, when a user adds a new landmark into a slice with aspline, the spline is immediately modified using precisely the samealgorithm of FIG. 39, namely operations 774 b-774 e. When a user moves alandmark, the spline is updated during the motion using operations 774b-774 e of the algorithm of FIG. 39. Adding landmarks (operations 774g-774 k of the algorithm of FIG. 39) is avoided during the motion phaseas it may lead to multiple updates during motions, resulting in too manypoints.

Once the contour lines or splines are successfully segmented from eachtarget image slice, the contour lines or splines are compiled asdiscussed above into a 3D mesh that may be used as an arthritic bonemodel 36 (see FIG. 1D) or bone models 22 (see FIG. 1C).

In one embodiment of the registration process discussed above, anoptimizer may be used during the registration process to maximizesimilarity between the open golden mesh and landmarks in the targetimage (and possibly edges image) by adjusting the parameters of a giventransformation model to adjust the location of reference imagecoordinates in the target image. In one embodiment, the optimizer for aregistration operation may use the transformed golden femur data fromthe previous registration operation as its initial approximation. Then,local optimization techniques may be used to search for a local optimumnear the initial starting approximation. This may be done so that anypotential matches farther away from the feature of interest (e.g., thefemur or tibia in a knee joint) reliably found in an earlier operationmay be eliminated.

In operation 770 a of FIG. 37, when optimizing the translationtransformation, a regular step gradient descent optimizer may be used byone embodiment. Other embodiments may use different optimizationtechniques.

To find a local minimum, parameter steps may be taken in the directionof the negative of the metric gradient (or the approximate gradient)over the transform parameter space at the current point. This generallyoptimizes the metric which typically has a local minimum when featuresof the reference image mapped into corresponding features of the targetimage have minimal misalignment).

In one embodiment, initial gradient step of 3 millimeters may bespecified, a relaxation factor may be set to 0.95 and a maximum of 50iterations may be used in the regular step gradient descent optimizationmethod to determine the parameters of the translation transformationthat results in minimal misalignment between the reference Open GoldenFemur mesh and the Landmarks in the target MRI scan.

In operation 770 b of FIG. 37, when optimizing the similaritytransformation, a regular step gradient descent optimizer may be usedagain by one embodiment. When applying the regular step gradient descentoptimizer to similarity transformation, the result and the convergencerate depend on the proper choice of parameters representing thesimilarity transforms. A good choice of parameters when used withgradient computations is such that a variation of every parameter by oneunit leads to approximately equal displacement of object points. Inorder to ensure similar displacement of points with respect to threerotational parameters in the similarity transform, the initial center ofrotation for the similarity transformation may be specified as thecenter of a bounding box (or minimum sized cuboid with sides parallel tothe coordinate planes) that encloses the feature (e.g., a bone)registered in the translation registration (e.g., operation 770 a inFIG. 37). For knee segmentation, scaling coefficients of approximately40-millimeters may be used for the scaling parameters when bringing therotational angle parameters together with translation parameters. Ascaling coefficient of approximately 40-millimeters may be used becauseit is approximately half the size of the bone (in the anterior/posteriorand medial/lateral directions) of interest and results in a point beingmoved approximately 40-millimeters when performing a rotation of oneradian angle. By the same reason a scaling coefficient of 40 millimetersmay be used in the similarity transform scaling parameter together withits translational parameters.

In one embodiment, an initial gradient step of 1.5 millimeters may bespecified, a relaxation factor may be set to 0.95 and a maximum of 50iterations may be performed in the regular step gradient descentoptimization method to determine the parameters of the similaritytransformation that results in minimal misalignment between thereference open golden template mesh and landmarks in the target MRIscan.

In operation 770 c of FIG. 37, when optimizing the affinetransformation, a regular step gradient optimizer may be used again byone embodiment. For knee bones, scaling coefficients of approximately 40millimeters may be used for the matrix coefficients variations whenbringing them together with translation parameters. An initial gradientstep of 1 millimeter may be specified, the relaxation factor may be setto 0.95 and a maximum of 50 iterations may be performed to determine theparameters of the affine transformation that results in minimalmisalignment.

In operation 770 e of FIG. 37, when optimizing the B-splinetransformation, a modified regular step gradient descent optimizer maybe used by one embodiment when searching for the best B-splinedeformable transformation. Namely, a combination of regular stepgradient descent optimizer with by coordinate descent may be used here.Rather than computing one gradient vector for the transform space andtaking a step along it, a separate gradient may be computed for everyB-spline transform node. In one embodiment, order three B-splines (withJ×K×L control nodes) may be used and J×K×L gradients may be computed,one for each control point. At every iteration, each of the spline nodesmay be moved along its respective gradient. This may enable fasterconvergence of the optimization scheme. A relaxation factor of 0.95 maybe used for each spline node. A an initial gradient step ofone-millimeter may be set for every B-spline grid node, and a maximum of50 iterations may be used in the regular step gradient descentoptimization method to find the parameters of the B-splinetransformation that provides minimal misalignment of the open goldenfemur mesh and landmarks and feature edges in the target MRI scan.

FIG. 23 depicts a flowchart illustrating one method for generatingspline curves outlining the surface of an object of interest in eachtarget MRI slice (e.g., as discussed above with respect to operation 772of FIG. 36) after the transformed golden femur mesh is found inoperation 770 e in FIG. 37. Initially, operation 470 intersects thetransformed golden femur mesh with a slice of the target scan data. Theintersection defines a polyline curve of the surface of the feature(e.g., bone) in each slice. Two or more polyline curves may be generatedin a slice when the bone is not very straightly positioned with respectto the slice direction.

A polyline curve is a piecewise linear approximation to a curved featureshape. Generally, this curve should be easy to manipulate with a set ofcontrol points. The polyline curve may have many segments, making itmore difficult to manipulate the polyline curve (e.g., during operation254 or 260 of FIG. 6). One embodiment may generate one or more Kochaneksplines from the polyline curve. Each spline typically has a smallernumber of control points and typically fits the polyline curve withabout 0.3-millimeter deviation. See section b. of this DetailedDescription for a detailed discussion regarding spline generation.

As discussed above, in one embodiment, the output of the segmentationmay be a triangular mesh (e.g., a 3D surface model) of the segmentedbone(s) of a joint (e.g., the femur and tibia of a knee joint). The meshgenerated generally represents a watertight surface that closely followsthe segmentation contour curves of the slices, smoothly interpolatesbetween the segmentation contour curves, and may have a low triangularcount. See section b. of this Detailed Description for a detaileddiscussion regarding mesh generation and the manual adjustment ofsegmentation splines.

The 3D surface models of the lower end of the femur and the upper end ofthe tibia of a patient's knee may be used to create arthroplasty jigsand/or implants. For example, the models may be used to create femur andtibia jigs that can be used with a patient's femur and tibia asdisclosed in the various U.S. Patent Applications incorporated byreference herein in this Detailed Description and filed by Park and Parket al. The automatic or semi-automatic processes described herein forsegmentation of image data to generate 3D bone models may reduce theoverall time required to perform a reconstructive surgery to repair adysfunctional joint and may also provide improved patient outcomes.

Although the present invention has been described with respect toparticular embodiments, it should be understood that changes to thedescribed embodiments and/or methods may be made yet still embraced byalternative embodiments of the invention. For example, certainembodiments may operate in conjunction with a MRI or a CT medicalimaging system. Yet other embodiments may omit or add operations to themethods and processes disclosed herein. Accordingly, the proper scope ofthe present invention is defined by the claims herein.

What is claimed is:
 1. A method for pre-operative planning for anarthroplasty procedure for a patient bone in an arthroplasty targetarea, the method comprising: receiving scan data associated with thepatient bone corresponding to the arthroplasty target area, the scandata captured using a medical scanning system; generating a model of thepatient bone from a golden template that has been modified using thescan data, the golden template being associated with an exemplary bone;generating preoperative planning data from the model, the preoperativeplanning data corresponding to an arthroplasty implant to repair atleast a portion of a joint associated with the patient bone; andgenerating a treatment based on the preoperative planning data, thetreatment including a bone resectioning plan for the arthroplasty targetarea to prepare the arthroplasty target area based on the arthroplastyimplant.
 2. The method of claim 1, wherein the scan data includes imageslices of the patient bone.
 3. The method of claim 1, wherein the modelof the patient bone is generated based on landmarks associated with thepatient bone in the scan data.
 4. The method of claim 1, wherein theexemplary bone is a bone type that is the same as the patient bone andis representative of what is considered to be non-diseased andnon-damaged with respect to the bone type.
 5. The method of claim 1,wherein the golden template is modified by mapping the golden templateto the scan data.
 6. The method of claim 1, wherein the model ismodified to represent a non-diseased state of the patient bone.
 7. Themethod of claim 1, wherein the bone resectioning plan includes at leastone of cutting, drilling, reaming, or resurfacing the arthroplastytarget area.
 8. The method of claim 1, wherein one or more machines areguided with a computer based on the bone resectioning plan.
 9. A methodfor pre-operative planning for an arthroplasty procedure for a patientbone in an arthroplasty target area, the method comprising: receivingscan data associated with the patient bone corresponding to thearthroplasty target area, the scan data captured using a medicalscanning system; identifying one or more landmarks in the scan data;obtaining golden bone data associated with an exemplary bone; generatinga patient model of at least a portion of the patient bone by modifyingthe golden bone data using the one or more landmarks in the scan data;generating a non-diseased state of the patient model; and generatingpreoperative planning data from the non-diseased state of the patientmodel, the preoperative planning data corresponding to an arthroplastyimplant to repair at least a portion of a joint associated with thepatient bone.
 10. The method of claim 9, wherein the repair includesrestoring the joint to at least one of natural alignment or zeromechanical axis alignment.
 11. The method of claim 9, wherein the scandata includes image slices of the patient bone.
 12. The method of claim9, wherein the golden bone data includes a model of the exemplary bone.13. The method of claim 9, wherein the golden bone data is modified bymapping the golden bone data to the scan data using the one or morelandmarks.
 14. The method of claim 9, wherein bone resection planningdata is generated based on the preoperative planning data and specifyinga treatment to prepare the arthroplasty target area for the arthroplastyimplant.
 15. A method for pre-operative planning for an arthroplastyprocedure for a patient bone in an arthroplasty target area, the methodcomprising: receiving scan data associated with the patient bonecorresponding to the arthroplasty target area, the scan data capturedusing a medical scanning system; generating a patient model of at leasta portion of the patient bone from golden bone data that has beenmodified by mapping the golden bone data to the scan data using one ormore landmarks, the golden bone data being associated with an exemplarybone, the patient model corresponding to a non-diseased state of thepatient bone; and generating preoperative planning data from the patientmodel, the preoperative planning data corresponding to an arthroplastyimplant to repair at least a portion of a joint associated with thepatient bone.
 16. The method of claim 15, wherein the golden bone dataincludes a model of the exemplary bone.
 17. The method of claim 15,wherein the patient model is generated using one or more curves modifiedto match the one or more landmarks.
 18. The method of claim 17, whereinthe one or more curves are generated from a golden bone mesh of thegolden bone data, the golden bone mesh being deformed to match the oneor more landmarks.
 19. The method of claim 15, wherein bone resectionplanning data is generated based on the preoperative planning data andspecifying a treatment to prepare the arthroplasty target area for thearthroplasty implant.
 20. The method of claim 19, wherein one or moremachines are guided with a computer based on the bone resectioningplanning data.